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ABSTRACT 


In this thesis, the finite-element method was developed 
to numerically analyze heat transfer by laminar natural 
convection within a rectangular cavity, a classical fluid 
flow problem. A second auxiliary case study involving 
Couette flow was included to test the flexibility of this 
analysis technique. 

Analyzing heat flows experimentally was also explored 
utilizing holographic interferometry. Specific problems 
encountered during this phase of research are presented 


with appropriate comments. 


TABLE OF CONTENTS 


NOMENCLATURE ---------------------------------------- 7 
I. INTRODUCTION ---------------------------------- 9 
A. HOLOGRAPHIC INTERFEROMETRY ---------------- 9 
B. CONCEPT AND HISTORY OF THE FINITE- 
ELEMENT METHOD --------------0--- seen n---- 12 
II. FUNDAMENTAL THEORY OF FINITE-ELEMENT METHOD --- 17 
A. VARIATIONAL PRINCIPLE --------------------- 23 
B. FINITE-ELEMENT METHOD AND THE RITZ 
TECHNIQUE --9rernre nme rte nnn rec rerne rent e- 24 
C. METHOD OF WEIGHTED RESIDUALS 
(GALERKIN'S METHOD) ----------------------- 29 
IIIT. ANALYSIS OF CONVECTIVE HEAT TRANSFER BETWEEN 
PARALLEL PLATES ---erese re renner nna rere nn 37 
a weer Cel REOULLO =———-————————————— == 41 
A. FINITE-ELEMENT ANALYSIS OF CONTINUUM ------ 41 
1. Discretization of the Continuum ------- 41 
2. Selection of Interpolation Functions -- 43 


3. Determination of the Elemental 


Properties -----en cnc ce rrr rete n nnn n-- 45 

B. DERIVATION OF ELEMENT MATRICES ------------ 51 
C. STRUCTURE OF COMPUTER PROGRAMS FOR 

FLOW ANALYSIS ------------- renner rer rrr nn- 62 


D. NUMERICAL RESULTS ------------------------- 71 


DeEEEPERTMENTAL PROCEDURE --~-----------------=me-—- iG 


A. ARRANGEMENT OF TEST APPARATUS ------------- 74 

B. HOLOGRAPHIC INTERFEROMETRY APPLICATIONS --- 80 

VI. CONCLUSIONS ------------------ portent enn ------- 84 
APPENDIX A FIGURES -------------------------------- 86 


APPENDIX B BRIEF REVIEW ON CALCULUS OF VARIATIONS - 98 


COMPUTER OUTPUT ween n men rrr re rn rn ret errr n= 102 
eee mee FLOW PROG LEM ==---—--—-————— 102 
B. STEADY STATE HEAT TRANSFER PROBLEM ---------- 108 
GOPIPUTER PROGRAMS 9-9-3 99033 9 = === Ss LLY, 
ae, GOUlM Tin LOW PROBEEM ~-<-=—<-——--=<---—-=->-_-—— ey, 
B. STEADY STATE FLUID MECHANICS PROBLEM -------- 134 
C. TIME-DEPENDENT FLUID MECHANICS PROBLEM ------ 153 
POP GENER ENG RS - c= —=—=---.2--~=-Ss"-s—2- 7 oa aoe T9353 
PEMA DEOLREBULION LEST ----="<-----==--=~~ oes os 195 


Cp 


NOMENCLATURE 


specific heat of fluid 
width of the enclosure 
area of element triangle 


gravitational acceleration 
gBL>(Ty-Te) 
VY 


Grashof number in L direction = 
height of the enclosure 
natural cooveiiates 
interpolation functions 
pressure (either wall or fluid) 
Prandtl number - £- 


3 
Rayleigh number in D direction = gBD (THT) 
a 


3 
Rayleigh number in L direction = 8BL (Ty-T¢) 
yr 


temperature (either wall or fluid) 
time relative to beginning of solution 
temperature at cold wall 

temperature at hot wall 


mean temperature of fluid = (TytTo) 
2 


velocity in x-direction 


‘velocity in y-direction 


Pd NRX BW R= * 


independent coordinate in horizontal direction 
independent coordinate in vertical direction 


thermal diffusivity of fluid = —% 
PP Cp 


coefficient of thermal expansion of fluid 
thermal conductivity of fluid 

dynamic viscosity of fluid 

kinematic viscosity of fluid = <A 
fluid density 


domain of integration for element (e) 


I. INTRODUCTION 


A. HOLOGRAPHIC INTERFEROMETRY 

The phenomenon of interference has had a considerable 
influence on the development of physics. Thomas Young's 
observation and explanation of the interference of the 
beams through two holes provided the basis for Fresnel's 
wave theory of light and the same experiment has been used 
as the foundation of modern coherence theory. 

Derived from interference is the technique of inter- 
ferometry, now one of the important methods of experimental 
physics. The father of visible-light interferometry was 
A. A. Michelson, who was awarded in 1907 the Nobel prize in 
physics for "his optical instruments of precision and the 
spectroscopic and metrological investigations he has executed 
with them.'' Applications to other spectral regions were 
more recent: the first use of interferometry in radio 
astronomy was reported in 1947, and infra-red interference 
spectroscopy was successfully employed some thirteen years 
later. 

Ever since the wave ees of light was generally 
accepted, interferometry has been the primary method for 


making measurements with great accuracy. The very small 


wavelength of light, on the order of 5x107-cm, and the fact 
that interferometric means are available for detecting 
changes of only a small fraction of this length, indicates 
the degree of accuracy which can be achieved. The wide- 
spread applications of the method attest to its general 
| usefulness. Interferometry is used for testing optical 
components, optical gauging of machine tools, studying air 
flow in wind tunnels, and standardizing the fundamental 
units of length. Therefore it is understandable that any 
fundamental improvement or innovation in this interferometric 
technique would find many applications over a wide field. 
Holographic interferometry is just such an innovation. 
Holography may be described as a photographic technique in 
which the amplitude and phase characteristics emanating 
from a coherent light source are recorded and later repro- 
duced. This reproduction assumes the form of a three- 
dimensional image of the original subject. Holography has 
widened the scope of interferometry to such a degree that 
holographic interferometry is now considered a standard tool 
in engineering laboratories all over the world. 
Conventional interferometry can be utilized to make 
measurements on highly polished surfaces of relatively 


simple shape. Holographic interferometry extends this 
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range by allowing measurements to be made on three- 
dimensional surfaces or arbitrary shape and condition. 
A roughly processed machine part can now be measured to 
optical tolerance. Furthermore, with the holographic 
technique a complex object can be examined interferometri- 
cally from many different perspectives, because of the 
three-dimensional nature of the hologram. A single inter- 
ferometric hologram is equivalent to many observations 
with a conventional interferometer. This property is 
especially useful for observations of such things as fluid 
flow in a wind tunnel. A third departure of holographic 
interferometry from conventional interferometry is that an 
object can be interferometrically examined at two different 
times; one can detect with wavelength accuracy any changes 
undergone by an object over a period of time. The present 
object can thus be compared with itself as it was at an 
earlier time. This is a great advantage in many fields. 
For example, a large lens can be tested before and after 
mounting. Similarly, with the aid of pulsed lasers, a 
machine part can be interferometrically compared with it- 
self statically as well as dynamically. 

Methods of holographic interferometry include single- 


and double-exposure as well as pulsed laser interferometry. 
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In this thesis, only single-exposure holographic interfer- 
ometry was considered since it corresponds to real-time 
. interferometry, that is, a method which allows one to 


observe changes in a subject as they actually occur. 


B. CONCEPT AND HISTORY OF THE FINITE-ELEMENT METHOD 

One must often resort to numerical procedures in order 
to obtain quantitative approximate sealants to linear and 
nonlinear problems in continuum mechanics. However, 
regardless of the initial assumptions and the methods 
used to formulate a problem, if numerical methods are 
employed in evaluating the results, the continuum is, in 
effect, approximated by a discrete model in the solution 
process. This observation suggests a logical alternative 
to the classical approach, namely, represent the continuum 
by a discrete model at the onset. One such approach, 
based on the idea of piecewise approximating continuous 
fields, is referred to as the finite-element method. Its 
simplicity and generality make it an attractive candidate 
for applications to a wide range of engineering problems. 

Classically, the analysis of continuous systems often 
began with investigations of the properties of small differ- 
ential elements of the continuum under investigation. 


Relationships were established among mean values of various 
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quantities associated with the infinitesimal elements, and 
partial differential equations or integral equations 
governing the behavior of the entire domain were obtained 
by allowing the dimensions of the elements to approach 
zero as the number of elements became infinitely large. 

In contrast to this classical approach, the finite- 
element method begins with investigations of the properties 
of elements of finite dimensions. The equations describing 
the continuum may be employed in order to arrive at the 
properties of these elements, but the dimensions of the 
elements remain finite in the analysis, integrations are » 
replaced by finite summations, and the partial differential 
equations of the continuous media are replaced, for example, 
by systems of algebraic or ordinary differential equations. 
The continuum with infinitely many degrees of freedom is 
thus represented by a discrete model possessing a finite 
number of degrees of freedom. Moreover, if certain com- 
pleteness conditions are satisfied, then, as the number of 
finite elements is increased and their dimensions are 
decreased, the behavior of the discrete system converges 
to that of the continuous system. A significant feature 
of this procedure is that, in principle, it is applicable 


to the analysis of finite deformations of materially 
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nonlinear, nonhomogeneous bodies of any geometrical shape 
with arbitrary boundary conditions. 

The practice of representing a structural system by a 
collection of discrete elements dates back to the early 
days of aircraft structural analysis, when wings and fuse- 
lages, for example, were treated as aseenblages of stringers, 
skins, and shear panels. By representing a plane elastic 
solid as a collection of discrete elements composed of bars 
and beams, Hennikoff /1941/ introduced his "framework 
method," a forerunner to the development of general dis- 
crete methods of structural mechanics. Topological prop- 
erties of certain types of discrete systems were examined 
by Kron /1939/, who developed systematic procedures for 
analyzing complex electrical networks and structural 
systems. Courant /1943/ presented an approximate solution 
to the St. Venant torsion problem in which he approximated 
the warping function linearly in each of an assemblage of 
triangular elements and proceeded to formulate the problem 
using the principle of minimum potential energy. Courant's 
piecewise application of the Ritz method involves all the 
basic concepts of the procedure now known as the finite- 
element method. In 1954, Argyris and his collaborators 


began a series of papers in which they developed certain 
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generalizations of the linear theory of structures and 
presented procedures for analyzing complicated discrete 
structural configurations in forms easily adapted to the 
digital computer. 

The formal presentation of the finite-element method 
together with the direct stiffness method for assembling 
elements was attributed to Turner, Clough, Martin, and 
Topp £19567, who employed the equations of classical 
elasticity to obtain properties of a triangular element 
for use in the analysis of plane agPese problems. It was 
Clough /19607, who first used the term "finite elements" 
in a later paper devoted to plane elasticity problems. 

Concepts of the method became more understandable 
after 1963 when Besseling /1969/7, Melosh /1970/, Fraeys 
de Veubeke /1971/, and Jones /1972/ recognized that the 
finite-element method was a form of the Ritz technique and 
demonstrated its generality for handling elastic continuum 
problems. In 1965, the finite-element method received an 
even broader interpretation when Zienkiewicz and Cheung 
[1973/7 reported that it was applicable to all field prob- 
lems which could be cast into variational form. During 
the late 1960's and early 1970's, while mathematicians 


were working on eStablishing errors, bounds, and convergence 
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criteria for finite-element approximations, engineers and 
other appliers of this same method were also studying 
similar concepts for various problems in the area of solid 
mechanics.: 

Although a major portion of the literature written to 
date on the finite-element method deals with static and 
dynamic structural analysis, there has been a continuing 
steady increase in the number of applications in other 
fields. The goal of this thesis was to develop a computer 
program, utilizing the finite-element method, which could 
accurately analyze laminar natural convection within a 
vertical rectangular enclosure. The program should be able 
to properly analyze axisymmetric as well as two-dimensional 


flows. 
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II. FUNDAMENTAL THEORY OF FINITE-ELEMENT ANALYSIS 


In this section the fundamental theory on which the 
thesis was based is presented. Highlighted topics include 
the variational principle, some basic concepts of finite- 
element analysis and the Ritz technique, and finally the 
method of weighted residuals featuring the Galerkin crite- 
rion. The variational principle and the Galerkin method 
are looked at in detail in regards to the derivation of 
finite-element equations. 

The finite-element method envisions a solution region 
as built up of many small, interconnected subregions or 
elements. Such a model of a problem gives a piecewise 
approximation to the governing equations. The basic premise 
of the finite-element method is that a solution region can 
be analytically modeled or approximated by replacing it 
with an assemblage of discrete elements. These finite- 
element discretization procedures reduce the problem to 
one of a finite number of unknowns by dividing the solution 
region into elements and by expressing the unknown field 
variable in terms of assumed approximating or interpolation 


functions within each element. The interpolation functions 


iy] 


are defined in terms of the value of the field variables 

at specified points called nodes or nodal points. Nodes 
usually lie on the element boundaries where adjacent elements 
are considered to be connected. In addition to boundary 
nodes, an element may also have a few interior nodes (al- 
though this was not the case in the choice of linear and 
quadratic trinagular elements utilized in this thesis). 

The nodal values of the field variable and the interpolation 
functions for the elements completely define the behavior 

of the field variable within the elements. For the finite- 
element representation of a particular problem, the nodal 
values of the field variable become the new unknowns. Once 
these unknowns are found, the chosen interpolation functions 
define the field variable throughout the assemblage of 
elements. 

Clearly, the nature of the solution and the degree of 
approximation depend not only on the size and number of the 
elements used, but also on the interpolation functions 
selected. As one would expect, functions cannot be arbi- 
trarily chosen since certain compatibility conditions must 
be satisfied. Often such functions are selected so that 
the field variable and/or its derivatives are continuous 


across adjoining element boundaries. Another important 
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feature of the finite-element method which sets it apart 
from other approximate numerical methods is its ability to 
formulate solutions for individual elements before putting 
them together to represent the entire problen. 

The finite-element method has gained much popularity 
and has been utilized extensively in recent years because 
it has, in general, several outstanding advantages. These 
are the following: 

1. Non-homogeneous configurations may be treated with 
relative simplicity. 

2. The elements can be graded in size and shape to 
follow boundaries of arbitrary shape. 

3. Once a computer program has been developed, problems 
of the same variety can be solved simply by supplying the 
computer with appropriate data. 

There are at least three distinct approaches one may 
employ in order to obtain finite-element equations of a 
particular system. In order of increasing versatility they 
are: (1) the direct approach, (2) the variational principle, 
and (3) the weighted residuals approach. 

The direct approach can be used only for relatively 
simple problems in which discrete elements may be easily 


identified. Once these elements have been selected, direct 
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physical reasoning is introduced to establish the element 
equations in terms of pertinent variables. The final step 
is then to combine the element equations to form the ge¥een- 
ing equations of the complete system. 

A detailed explanation of the remaining two approaches 
will be given in Subsections A, B and C to follow. 

Whichever one of these three particular approaches is 
utilized, the finite-element method follows a systematic 
step-by-step process when applied to continuum problems. 
They are: 

1. Discretize the continuum. 

The entire flow region under study is divided into 
a series of subregions or elements assumed to be inter- 
connected at a finite number of nodal points; thus a program 
originally exhibiting an infinite number of degrees of free- 
dom is made finite. The elements used can be triangular, 
rectangular, or almost any shape. Also, information must 
be fed into a computer giving global coordinates of the 
nodes and topology of the system. 
Finally, selection of which field variables are to 

be used to satisfactorily describe solution domain must 


be indicated at this point in the process. 
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2. Select interpolation functions, n, 6@) 

From the nodal values one represents the value of 
the field variable over the element by means of interpola- 
tion functions. Often, although not always, polynomials 
are selected as these functions because they are easy to 
integrate and differentiate. The number of nodes and the 
order of the interpolation polynomials are interrelated. 
The field variable itself may be a scalar, a vector, or a 
higher-order tensor. 

3. Find the element properties. 

Essentially, the problem is solved at the element 
level. The matrix equations expressing the properties of 
the individual elements are determined. This can be 
accomplished by any one of the three approaches previously 
mentioned: the direct method, the variational principle, 
or the weighted residual method. The approach used depends 
entirely on the nature of the particular problem. 

4. Assemble the element properties to obtain the system 
equations. 

In this step, one combines the matrix equations ex- 
pressing the behavior of the elements to form the matrix 
equations expressing the behavior of the entire solution 


region or system. The matrix equations for the system 
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exhibit the same form as the equations for an individual 
element except that they contain many more terms because 
they include all the nodes. The basis for this assembly 
procedure stems from the fact that, at a node where elements 
are interconnected, the value of the field variable is the 
same for each element sharing that node. 

5. Solve the system equations. 

From the previous step, a set of simultaneous equa- 
tions are derived which can now be solved to obtain the 
unknown nodal values of the field variable. If these 
equations are linear, a number of standard solution tech- 
niques may be employed; if the equations are nonlinear, 
their solution is more difficult to obtain, but several 
alternative approaches do exist that lead to satisfactory 
results. 

6. Make additional computations if desired. 

It may be desired to use the solution of the system 
equations to calculate other important parameters, i.e., 
from the nodal values of the pressure, one might wish to 
calculate velocity distributions. 

It is worth making mention of the fact that several of 
the steps in the above process are essentially the same 


regardless of the type of problem (this thesis was devoted 
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to the fluid mechanics problem). Thus, only steps three 

(3) and six (6) might differ for any given situation, in 
that the equations describing the elements could vary. The 
other steps would be the same. This generality of the 
finite-element method is, without doubt, one of its greatest 


strengths. 


A. VARIATIONAL PRINCIPLE 

Often, continuum problems have different but yet equiva- 
lent formulations, such as a differential formulation and 
a variational formulation. In the differential case, the 
problem is to integrate a differential equation or a system 
of differential equations subject to given boundary condi- 
tions. In the classical variational formulation, the 
problem is.to find the unknown function or functions which 
extremize (maximize, minimize) or make stationary a func- 
tional or system of functionals subject to the same specified 
boundary conditions. The two problem formulations are 
equivalent because the functions that satisfy the differential 
equations and their boundary conditions also extremize or 
make stationary the functionals. This equivalence is 
apparent from the calculus of variations, which shows that 


the functionals are extremized or made stationary only when 
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one or more Euler equations and their boundary conditions 
are satisfied. Consequently, these equations are precisely 
the governing differential equations of the problem. To 
illustrate this duality concept, Appendix B provides a brief 
review and introduction to some basic ideas of the calculus 


of variations. 


B. FINITE-ELEMENT METHOD AND THE RITZ TECHNIQUE 

The Ritz technique is basically a procedure for trans- 
forming a continuous medium into an approximated lumped 
parameter system. A more qualitative definition would be 
that the Ritz method consists of assuming the form of the 
unknown adjustable parameters. From this family of trial 
or coordinate functions, that particular function which 
renders the functional stationary is then selected. The 
procedure is to substitute the trial functions into the 
functional and thereby express the functional in terms of 
the adjustable parameters. This functional is then differ- 
entiated with respect to each parameter, and the resulting 
equation is set equal to zero. If there are n_ unknown 
parameters, there will be n_ simultaneous equations to be 
solved for these parameters. By this means, the approximate 


solution is chosen from the family of assumed solutions. 
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This procedure does nothing more than give one the 
"best" solution from the family of assumed solutions. 
Clearly, then, the accuracy of the approximate solution 
depends on the choice of trial functions. These trial 
functions are required to be defined over the whole solu- 
tion domain and must satisfy at least some and usually all 
of the boundary conditions. Sometimes, if the general 
nature of the desired solution is known, the approximation 
can be improved by choosing the trial functions to reflect 
this nature. If, by chance, the exact solution is contained 
in the family of trial solutions, then the Ritz technique 
gives the exact solution as expected. Generally, the 
approximation improves as the size of the family of trial 
functions and the number of adjustable parameters increase. 
If the trial functions are part of an infinite set of 
functions that are capable of representing the unknown 
function to any degree of accuracy, the process of including 
more and more trial functions leads to a series of approxi- 
mate solutions which converge to the true solution. Often 
a family of trial functions is constructed from polynomials 
of successively increasing degree. 

To illustrate the Ritz technique, consider the following 


simple example. Suppose it is desired to find the general 


ZD 


function @(x) satisfying 
d*¢ = -£(x) 
dx2 


with boundary conditions of @(a)=A and @(b)=B specified. 
It is assumed that f(x) is a continuous function in the 
closed interval /a,b/. This problem is equivalent to finding 


the function $(x) that minimizes the functional 


a 
1(9) = 54 | aah? - £060 | dx 
b 


x 
2 
SEs of the fox 1(¢) = S F(x,6, 0,, £(x))dx 
sia | 
Ignoring the fact that this problem possesses an exact 
solution, we will attempt to find an approximate solution. 
According to the Ritz method, the desired solution can be 


assumed to be approximately represented in fa,b/ by a con- 


bination of selected trial functions of the form 
Q(x) = Cy Y (x) +C, Fo(x)t... AC, 7 (x) , aSxSb 


where the n_ constants C. are the adjustable parameters to 
be determined. The trial functions should be selected so 
that the expression for @(x) satisfies the boundary condi- 


tions regardless of the choice of the constants C,. Using 
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polynomials is a simple and convenient way of constructing 


the trial functions. Therefore 


B(x) a7 (x-a) (x-b) (CC xtC gx +... HC a, 
n 


is a possible series of trial functions. When this approxi- 
mate expression for @(x) is substituted into the functional 
to be minimized, and after the integration has been carried 


out, the functional is of the form 


Since the Cc. are required to be chosen such that they mini- 
mize I, employing differential calculus, the following 


partial differential equations are formulated 


oI 19,91 -0,.,.., 91 =9 
OC} a Co dC, 


These n equations are then solved for the n parameters 
C.; and the accuracy of the approximate solution depends on 
the number of C's used in the trial function. Generally, 
as n increases the accuracy improves. To assess the 
improvement in accuracy as more C's are utilized, the prob- 
lem is solved repeatedly by taking successively more terms 


in the approximation, that is 
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X 


O,(=) % (x-a)(x-b) Cy 


Bo(x) AS (x-a)(x-b)(C,4C,x) 


ae 
$3(x) LD (x-a) (x-b) (C,4C,x4C,x”) 


and so on. By comparing the results at the end of each 
calculation, the effect on accuracy of adding more terms 
can be estimated. 

The finite-element method and the Ritz technique are 
essentially equivalent. Each method uses a set of trial 
functions as the starting point for obtaining an approxi- 
mate solution; both methods take linear combinations of 
these trial functions; and both methods seek the combination 
of trial functions that renders a given. functional stationary. 
The major difference between these approximating methods 
stems from the fact that the assumed trial functions in the 
finite-element method are not defined over the entire solu- 
tion domain, and they must satisfy not just any boundary 
conditions, but only certain continuity conditions and then 
only sometimes. Since the Ritz technique uses functions 
construed over the whole domain, it can be employed only 
for domains of relatively simple geometric shape. Also, 
these trial functions associated with the Ritz method are 


required to satisfy at least some and usually all of the 
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boundary conditions. In the finite-element method the 

same geometric limitations exist, but only for the elements. 
Due to the fact that elements with simple shapes can be 
assembled to represent exceedingly complex geometries, the 
finite-element method is a far more versatile tool than the 
Ritz technique. From a strict mathematical standpoint, the 
finite-element method is a special case of the Ritz technique 
only when the piecewise trial functions obey certain con- 
tinuity and completeness conditions that are stipulated 


over just the element alone. 


C. METHOD OF WEIGHTED RESIDUALS (GALERKIN'S METHOD) 

The third and final approach to the finite-element 
method involves a procedure that is more generalized and 
straightforward than either of its two predecessors. 

The relationship between the well-known Ritz technique 
and the finite-element method enables one to view the finite- 
element discretization procedure as simply another —— 
for finding approximate solutions to variational problems. 
In fact, these finite-element equations were shown to be 
derived by requiring that a given functional be stationary. 
This broad variational interpretation is the one most 


widely used to derive element equations, and it is the 


Zo 


most convenient approach whenever a classical variational 
statement exists for a given problen. 

However, applied scientists and engineers encounter 
practical problems for which classical variational principles 
are unknown. In these cases finite-element Hewhnaamecmetne 
still applicable, but more generalized procedures character- 
istic of the method of weighted residuals must be employed 
to derive the element equations. Through certain generaliza- 
tions, finite-element equations may be derived directly 
from the governing differential equations of the problem 
without reliance on any classical, quasi-variational, or 
restricted variational "principles." This procedure allows 
one to apply the finite-element method to almost all practical 
problems of mathematical physics. 

The method of weighted residuals is a technique for 
obtaining approximate solutions to linear and nonlinear 
partial differential equations. It offers still another 
means with which to formulate the finite-element equations. 
Applying the method of weighted residuals involves basically 
two steps. The first step is to assume the general func- 
tional behavior of the dependent field variable in some way 
so as to approximately satisfy the given differential equa- 


tion along with its associated boundary conditions. 
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Substitution of this approximation into the original differ- 
ential equation and boundary conditions then results in 
some error called a residual. This residual is required to 
vanish in some average sense over the entire solution domain. 
The second step entails solving the equation(s) resulting 
from step one and thereby specializing the general func- 
tional form to a particular function, which in turn becomes 
the approximate solution sought. 

To be more specific, the following typical problem is 
offered. Suppose it is desired to find an approximate 
functional representation for a general field variable @ 
governed by the differential equation 

(9) - £ = 0 (2.1) 

in the domain D bounded by the surface » : Se is a 
linear or nonlinear differential operator and the function 
f is a known function of the independent variables. Also, 
proper boundary conditions are assumed to be prescribed on 
y . The method of weighted residuals is now applied in 
two steps. First, the unknown exact solution @ is approxi- 
mated by Or eneee either the functional behavior of ¢ is 
completely specified in terms of unknown parameters, or 

the functional dependence on all but one of the independent 


variables is given while the functional dependence on the 
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remaining independent variable is left unspecified. Thus 


the dependent variable is approximated by 
ys m 
$xG=) NAC, (22) 
i=1 


where the N. are the assumed functions and the C. are either 
the unknown parameters or unknown functions of one of the 
independent variables. The m functions N; are usually 
chosen to satisfy the global boundary conditions of the 
system in question. When @ is substituted into equation 


2.1, it is unlikely that this equation will not be satisfied, 


that is, 
LG) -£ #0 
but in fact, 
£8) -£ =e 


where e is the residual or error that results from approxi- 
mating @ on ty The method of weighted residuals seeks to 
determine the m unknowns C. in such a way that the error 

e over the entire solution domain is small. This is 
accomplished by forming a weighted average of the error 

and specifying that this weighted average vanish over the 
solution domain. In other words, m linearly independent 


weighting functions, Wi; are chosen such that 


J (e6-e] W,dD= J et,a0-0, i=1,2,---,m (2.3) 


ay. 


The form of the error distribution principle expressed 
in equation 2.3 depends on the choice of weighting functions. 
Once these are specified, equation 2.3 represents of a set 
of m equations, which may be either algebraic or ordinary 
differential. The second step is to solve for the C,'s and 
hence obtain an approximate representation of the unknown 
general field variable @ via equation 2.2. There are many 
linear problems and even some nonlinear problems for which 
it can be shown that, as m~o, ig, but, in general, 
studies of convergence and error bounds are scarce. 

Due to the broad choice of weighting functions or error 
distribution principles than can be used, a variety of 
weighted residual techniques are likewise available. The 
error distribution principle most often utilized to derive 
finite-element equations in the field of aeronautics is 
known as the Galerkin criterion, or Galerkin's method. 
Here, the weighting functions are chosen to be the same 
as the approximating functions employed to represent @, 


that is, W,=N, for i=1,2,....,m. Therefore Galerkin's 


i 
method requires that 


([e@-<] N,dD = 0 (2.4) 
D 


In the preceding section pertaining to the Ritz technique, 
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it was assumed that the entire solution domain was being 
dealt with. However, because equation 2.1 holds for any 
point in this region, it also holds for any collection of 
points defining an arbitrary subdomain or element of the 
whole domain. Consequently, attention may be focused 
directly on an individual element by means of a local approxi- 
mation analogous to equation 2.2, but being defined as valid 
for only one element at a time. Now the finite-element 
representations of a general field variable become available. 
The ——s N; become what are known as the interpolation 
functions n, fe) defined over the wae and the Cc, are 

the undetermined parameters, which may be the nodal values 
of the field variable or its derivatives. Then, from 
Galerkin's method, the equations governing the behavior of 


an element of the solution domain may be written as 


S$ [ecole -20] n, (©) a6 = 0, i=1,2,....,r (2.5) 
D 


where, as before, the superscript (e) restricts the range to 
one element, and 

gfe) «a be wee) | fg} Ce) 

¢(e) = forcing function defined over element (e) 


r = number of unknown parameters assigned to the 
element. 
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There exists a set of equations similar to equation 2.5 
for each element of the whole assemblage. Prior to assembling 
the system equations from the individual element equations, 
it is required that the choice of approximating functions N. 
guarantee the interelement continuity along the boundary 
necessary for the assembly process. If the field variable 
is continuous at element interfaces, then C° continuity 
exists; if, in addition, first derivatives of the variable 
are continuous, cl continuity is said to occur; if second 
derivatives are also continuous, a region of c? continuity 
exists; and so on. This is the standard definition and 
notation utilized for expressing the degree of continuity 
of a field variable at element junctions. The higher the 
order of continuity required in the solution, the narrower 
one's choice of interpolation functions becomes. 

With the above definition of continuity in mind, the 
compatibility and completeness requirements for such inter- 
polation functions may be stated. If the functions appearing 
under the integrals in the element equations contain deriva- 
tives up to the (ntl)th order, then the following stipulations 
must be satisfied for assurance of convergence as the element 


Size decreases. 
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Compatibility requirement: At element interfaces, ‘Ce 
continuity must exist. 


Completeness requirement: Within an element, cntl co 


n- 
tinuity must exist. These requirements hold regardless of 
whether the element equations (integral expressions) were 
derived using the variational technique or the Galerkin 
methed. For this thesis, n was taken to have a value of 
zero. 

Integration by parts is a convenient way to introduce 
the natural boundary conditions that must be satisfied on 
some portion of the system exterior or boundary. Although 
the boundary terms containing these imposed conditions 
appear in the equations for each element, during the assembly 
of the element equations only the boundary elements give 
nonvanishing contributions. After the assembly process has 
been completed, the fixed boundary conditions (i.e., speci- 
fied velocity, pressure or temperature) are conveniently 


introduced to help simplify the final matrix form of the 


finite element equation. 
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III. ANALYSIS OF CONVECTIVE HEAT 
TRANSFER BETWEEN PARALLEL PLATES 
The transfer of heat energy across a fluid layer is 
accomplished, in general, through the mechanisms of con- 
duction, convection and radiation. This last phenomenon 
ite usually a function of the fluid enclosed between the 
surfaces and the nature, temperature and configuration of 
the enclosing boundaries. Radiation takes place independ- 
ently of the conduction and convection as long as there is 
no absorption by the fluid, and therefore under these con- 
ditions it can be omeenaaues separately. The phenomena of 
conduction and emul are closely interdependent and | 
are usually analyzed together. Buoyancy forces result from 
differences in density within the fluid and are caused by 
heat transfer to or from this fluid. Natural convection 
may then be thought of as fluid motion of the system due 
to the activation of these buoyance forces. In a two- 
ieiaititiiens plane, such heat transfer across a vertical, 
enclosed fluid layer is a function of the Grashof number, 
the Prandtl number and the fluid layer height-to-width 


ratio (L/D). 


3/ 


Natural convection plays a very important role in 
materials processing at high temperatures where agitation 
by other means is impracticable, or where the existence 
of temperature gradients is an inherent characteristic of 
the system. 

The steady convective motion of a lubricating fluid 
contained within a long, rectangular enclosure was investi- 
gated. Holographic interferometry and numerical approxima- 
tion were the experimental and theoretical analysis tools, 
respectively. | 

The two vertical walls of the enclosure were held at 
different temperatures, and the top and bottom were deemed 
perfect insulators (Figure 1). It was considered that the 
length of the enclosure (7 inches) was sufficiently long 
in the direction normal to the plane of Figure 1 for the 
motion to be assumed two-dimensional. Another assumption 
made was that the fluid motion was laminar. Experimental 
evidence indicates that such an assumption is valid provided 
the Rayleigh number based on cavity height is less than 
about 10° (Ra, in this study was calculated to be 1.018 x 
10’). Using this value and a value of the Prandtl number 
of 1.0755 x 10, determined from the ratio of kinematic 


viscosity to thermal diffusivity of the fluid, a system 
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Grashof number of 946.4 was calculated. The temperatures 
of the vertical walls x=0 and x=D were defined to be Ty 
and To respectively. If (T,-T¢) in degrees Fahrenheit is 
sufficiently small with respect to Te; the Boussinesq 
approximation may be introduced which neglects density 
variations in inertia terms of the equations of motion, 
but retains it in the buoyancy term. One final assumption 
was made that all other relevant thermodynamic and trans- 
port properties were independent of temperature and that 
compressibility and viscous dissipation effects were 
negligible. 

The problem now was to find the time and spatial 
dependence of the velocities and the temperatures within 
the system. 

The governing differential equations expressing con- 
servation of mass, momentum (both in x- and y-directions) 


and energy were 


SU + S% = 0 (3.1) 
du Ou Ou d2u , d *u OP 


+gB(T-T.,) | (3.2) 


° 2 
m4 Vv 1 
Si tuft + - SS +  -ST O 


OT 
Ae ay axz dy? 


i, vot = x (ST + 9°%) (3.4) 


The solution to the foregoing set of dynamic equations must 


satisfy the following boundary conditions on the walls, 


ue =V5,™ 0, no velocity on any of the four walls 
T=Ty or To given on the two vertical walls 


P=PAaTMOS > also given on the two vertical walls. 
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IV. THEORETICAL RESULTS 


A, FINITE-ELEMENT ANALYSIS OF THE CONTINUUM 
1. Discretization of the Continuum 

Since the fundamental premise of the finite-element 
method is that a continuum or solution domain of arbitrary 
shape can be accurately modeled by an assemblage of simple 
shapes, most finite elements are geometrically simple also. 
This statement especially pertains to the choice of the 
triangular-shaped element which would represent the unknown 
System parameters in this study, that is, the velocity, 
temperature and pressure. The main reason behind this choice 
was the fact that the three-node flat triangular element is 
the simplest two-dimensional element available, and hence 
an assemblage of triangles could always depict a two- 
dimensional domain with any number of straight sides. The 
solution domain in this problem was the vertical rectangular 
enclosure, a relatively simple-shaped continuum which posed 
no problem for the triangular elements. Twelve (12) ele- 
ments were utilized to represent the 8.5 inch by 1.875 inch 
area. They were interconnected to each other and the bound- 
ary at a total of thirty-five (35) nodal points, of which 


twelve (12) were corner nodes (Figure 2). 
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Arriving at this figure of thirty-five nodes was 
not an arbitrary process. Once pressure was chosen to be 
linearly approximated, system velocities and temperature 
were required to assume polynomial approximation of one 
degree higher, or quadratic, if the highest solution 
accuracy was to be achieved. For triangular elements, a 
complete n th-order polynomial requires %(n+1)(n+2) nodes 
for its specification. Therefore, a lst order, or linearly 
approximated, polynomial is associated with a three-node 
triangle; and a quadratic er relates to a six-node 
triangle. 

The three-node elements, with their nodes on the 
corners, may be thought of as being superimposed onto the 
six-node elements. Such elements contain, in addition to 
the corner nodes, nodes located at the midpoint of each 
side of the triangle. Twelve triangular elements of the 
six-node variety may be interconnected to form the solu- 
tion domain shown in Figure 2; this domain possessing 
exactly thirty-five nodes. 

Each element (6-node and 3-node) specifies uniquely 
a complete polynomial of the order necessary to give C° 
continuity, and hence satisfy the completeness and compati- 


bility requirements for elemental assemblage. 
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Next, the distinction between local and global 
node-numbering had to be made. Since each element in the 
triangular mesh had six nodes, the local nodes were identi- 
fied as such by starting in the upper left hand corner of 
each element and numbering counterclockwise around the 
element. The global node system is a method for uniting 
these independent elements along with their nodes into one 
distinct entity. Figure 3 summarizes the relation between 
local and global numbering for four (4) such elements. 

This figure defines the system topology or the connectivity 
of the system. 
2. Selection of the Interpolation Functions 

In the preceding subsection it was mentioned that 
linear approximation was used for values representing nodal 
pressures, while both velocity and temperature varied in a 
quadratic fashion within the elements. Such a relationship 
was based on the governing equations of the system, in 
which the highest order of partial differential equations 
involving pressure was one, while partial derivatives of 
u, v and T existed up to second order. Therefore, choosing 
linear pressures required the remaining three nodal para- 


meters or field variables to take on quadratic approximation. 
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The functions employed to represent the behavior 
of these field variables within an element are known as 
interpolation or approximating functions. Their order 
within an element depends on the number of degrees of 
freedom assigned to that element. In this study, two 
different polynomial series were selected as the first 
and second order interpolation functions. Associated with 
these series were coefficients made up of generalized 
coordinates, that is, independent parameters which speci- 
fied the magnitude of the prescribed distribution for each 
field variable (u, v, P, T). These polynomials were repre- 


sented as follows 
P(x,y) Sac, (ete, (xtc, (e)y (4,1) 
for the linear pressure terms, and 
(x,y) (ac, (40, (atc, (©) y40,, (2) 224 
6,6 xytc, (€) y4 (4.2) 


with @ being a generalized quadratic field variable (either 


1 § 


u, v, or T in this case, and the superscript ‘e standing 


for element. 


ope 


The next step in the process was to solve for the 


generalized coordinate in terms of the as yet unknown 


c fe) 
i 
field variables. This gave the desired interpolation, but 

the form of the resulting equations was not convenient. 


As a final step then, the equations were rearranged until 


they appeared as 


P(x,y) Ce) ans? (x,y)Py4Ny* (x,y) Pot 


Na (x,y)P,-|_w_| fe} (4.3) 
and 


(x,y) (2) aN P(x, ¥)B, AN, M(x, y)ByANG (x,y) Og 


p p ® _| WG 
N, (x,y) B,ANs (x,y) BSANG (xy)H6=|_W_]f 6) 
(4.4) 
u V 1 ee : 
where N. , N, ,» and N, were the specific interpolation 
functions in equation (4.4) for this study and were all 


T 


equal in forn, i.e., N, =N.Y=N. =N.. 
iD oes ae 2 | 


3. Determination of the Elemental Properties 
In this thesis, the Galerkin method was utilized 
to determine the element properties. This procedure 
applied at a general node i of an isolated element 


becomes, in view of equations 3.1-3.4, 


45 


(e) (e) 
u 5 P=) V “i 
fac — al 


i 
a Su, y2ul a Bple) sa (qle)- 
f w| Bag are) = 4 SE scr a 


led ul _(e) dul)_ 2 SJ 7 
u _—— V yan . os dxdy=0 (4.6) 


Sm fe $v (e) aes (e wei vie ap? (e) sve) 
Ge Oy Poy Ox 


wold aul? | dxdy=0 (4.7) 
2r(e), J2r(e), (eae) Azle) 
Tt T -us&/ OT" 
be eS ye 8 
ao) OTS) ax | dxdy=0 (4.8) 


where W. (x,y) and H, (x,y) are the weighting or interpolation 
functions, which were taken as 

WAN, and H,=N,=L, (natural coordinates). 
The inertia terms in equations 4.6, 4.7 and 4.8 considerably 
increased the degree of difficulty of this fluid £low prob- 
lem when compared to an incompressible viscous flow without 
inertia. This is because the above mentioned equations are 
nonlinear, thereby forcing an iterate procedure to be 


introduced and repeated until the ule) , vie} , and Tie) 
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(e) 


values converged to the previous ue wey and Toe) 

n es 
solutions. The subscript n runs from zero to some posi- 
tive number at which the field variable passes a convergence 
test. 

Integrating each term of equations 4.5-4.8 by parts, 
and making use of the approximations of equations 4.3 and 
4.4, the following results on an elemental level were 


obtained 


i dxdy) fv} =0 (4.9) 
6 gi Mal =p aon (0 
non 
eSB dxdy fp} -2B Lx, | fr} 
ae© 
- ((n,| N | dxdy){ 24) = at) + f ni 
a? iL : xa 34} om Aa 


(4.10) 


{ 4Qu ote Nt +i QMi SEN) axay { v } 
A 
-~ + dxd - fa, | w J axayy {dv 
A) SAL wt tar} LN es 
= Ff N.Y*ds (4.11) 
c 
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fox (Qs Opt 4 2Ni OLN) axay .: } 


ae, 
- (N, | N_| dxdy) 47} = fyzhds 
‘ pe {gr 1 


Cc (4,12) 
where N,X*ds, N,Y*ds: and N,Z*ds are simply lumped-sum con- 
tour integrals niu introduce the natural boundary condi- 
tions for u, v and T respectively. These integral values 
were labeled QX, QY, and QZ in the computer program. The 
last term on the left hand side of equations 4.10-4.12 
represents the transient nature of the system. 

Finally, the element matrix equations were written 
by inspection from equations 4.9-4.12 and were of the 


general form 


Ea (e) g co). De | ek b} C= fz] (e) (4.13) 


where the square matrices [ K ] Ce) AND fs.) (e) are known 
as stiffness matrices, nhs column vectors { d } (e) and tg Hee) 
are the nodal field variable and time derivative vectors, 
respectively. The column vector ‘ R } (e) signifies the 
resultant nodal force vector for the element. In the actual 


computer program, the following identities were used 
fe] =[™]» [xe] - Lo] {9} -€e} fo} -fep 
and ( R } = t Rus} 
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and the element matrix equations were 


es] | [0] el; fo} | | fay 


ol E z el fol | | ey 


In the above assemblage, the individual matrix notation 


utilized was 


[| -KG.5)= és (QM ONS + Nt ONS) axay 
| Ky | 7=K."(i,5)= J (Ni Nj") dxay 
[Ky] “Ks (4,5) >| (RE nj") axay 
| Ky | =Ko(i,5)= vf (Qui N,*) dxdy 
[ K, ] -K3(4,5)= Jf (245 N, )dxdy 
“aS 


| co | -c0(4,5)= J w,ajaxay 
ow ; | 


Also, gBI_ is the u velocity forcing function, and r 

and s are the number of nodes where velocity (or temper- 
ature) and pressure are interpolated at, respectively. In 
this study, r=6 and s=3, therefore the element matrices were 
21x21 and the element column vectors 21xl. 

Once the matrix equations were compiled or assembled 
on the element level, assembling these properties to obtain 
the system equations in matrix form also was a relatively 
simple operation for the digital computer. In essence, the 


large square matrices were derived by systematically adding 
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together the contributions of each individual element matrix, 
and inserting prescribed nodal variables or boundary condi- 
tions where applicable. As was brought out in the theoreti- 
cal section of this thesis, the final assembly now became 

a system of ordinary differential equations resembling the 


same format as equation 4.13, i.e. 


(Jo) - [a] {3} fe) wu) 


The problem solution was completed when these equations 
were solved for the nodal parameters { 9} , subject to 


the discretized initial conditions. 


B. DERIVATION OF ELEMENT MATRICES 

Derivation of various element matrices, referred to 
previously as simply area integrals over the solution 
domain, will be discussed in this section. The evaluation 
of each matrix will be in terms of natural coordinates, 
that is, weighting functions relating the coordinates of 
the end nodes to the coordinate of any interior point 
belonging to the element. The weighting functions are not 
independent of one another, since their sum must equal 


unity, i.e. 


ee wee 
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where n is the number of external nodes of the element. 
This expression can be interpreted to mean that one and 
only one coordinate is associated with node i, having a 
unit value there and a zero value at every other node. As 
was previously mentioned ss other sections, a general 
triangular shaped element, such as sketched SS oes 
employed. Then by equation 4.15 

LitL5tL3 = 1 
A cartesian coordinate system is used since the fluid flow. 
is assumed to be two-dimensional. Similar results could be 
derived using cylindrical coordinates for an axisymmetrically 


shaped element. The original Cartesian coordinates of a 
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point in the element can now be linearly related to the 


new natural coordinates by the equations 


x= Ly x, TL, x, thx, 


and 


(4.16) 


(4.17) 


Solving for the natural coordinates in terms of the Cartesian 


coordinates gives 


—" 55 +—(a, +b, xtc, y) (4.18a) 
= 1 
Lo(x.y) = ga lagtboxteoy) (4, 18b) 
and finally 
= i 
L3(x,y) we (agtbaxte,y) (4.18c) 
where 
1 x Yi 
2QA= {1 X, Yo = 2 (area of triangle 1-2-3) 
gilt X, 3 


A1=X2¥3-X3V9, by=Yg~¥3, ©1=X37X2 
49=X3Y17°* 193» 2=V3-Y1>, C2=X1-*3 


23=%1¥o-XoVy> 3=¥y-Yo> C3™XQ-*y 


ao 


The interpolation functions N, for the linear pressures in 


terms of natural coordinates are merely 


E P P 
Nj =1,,N,=4,,N, =4, 


but those interpolation functions that relate to the 


velocities and temperatures stemming from quadratic approxi- 


mation possess the form 


N, = 2L,°-L, 

Ny = 4L,L, 

N, = 2L, -L, | cia 
N, = 4L5L, 

N, = 2L,-L, 

Ng = HyL, 


Another way of envisioning L, (x,y) for the triangular 

element is to consider it a ratio of areas. Figure 4 

shows how the natural coordinates, often called area 
coordinates, are related to areas. In this figure, when 

the point (x Y,) is located on the boundary of the element, 
one of the area segments vanishes and hence the appropriate 
area coordinate along that particular boundary is identically 


zero. For example, if om) is on line 1-2, then 
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Le 2 0 since A, = 0 


A 


There is also a convenient analytical method. for integrating 
area coordinates over the area of a triangular element and 
involves the formula 


1 Ly Lb,” date) = 61a! ¥ 


L AO. 8 
3 (X+8+8 +2)! 


(e) 


2A 
A 
A summation of values derived using this formula is pre- 


sented in Table 4.] for (K+A+3B) < 4. 


“x 8B 8 
lt lL, Ly L, anf) = A 


Table 4.1 
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With this preliminary work finished, the actual deriva- 
tion of the element matrices may now begin. In all, five 
matrices will be completely evaluated while one matrix, 
| ky | » will have only two of its terms derived, due to 
the extensive amount of time and paper needed to evaluate 
[| in total. Beginning with this above-mentioned 


matrix as it appeared in Subsection A, 


2 ON ion ONi ON} 
K (i,j) = V airs os. SF dxdy (4.20) 


nok 


where rowan is the elemental area representing the solution 
domain, and i and j both vary from one to six. Since 
lx, | is an array multiplying the nodal variables of 
velocity and temperature, it must be correlated with the 
quadratic interpolation functions of equation 4.19. For 


the point (1,1), equation 4.20 becomes 


K,(1,1) = vfcgu gel r ant 3%) axay (4.21) 
) 


(e 
where aie 
21-3 el = (4L,-1) = — (4L, -1) 
and 
ON, oly _ 
. = oes (4L, -1)= a (4L, -1) 


substituting these values into equation 4.21 
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2 2 
- by (au. -1) 7421 _ (4. - 2 | aed 
Ae 


yA 
b 2 
= gar CT) (16L, =<8L_+1) dxdy 
4h 2 es 
OQ 


employing Table 4.1 for these three cases above in which 


(X+tA+5) = 2, 1, and 0 respectively; plus the relation- 


ship that S dxdy =~ A 
a. 


2...2 
K,(1,1) = A EES Ges ie WES 
4A 12 3 


or finally 


2 
- ¥(P1 +¢ 
K, (1,1) - 


. (e) ¢ . d . 
where 
Shag | yi aySb | . 
ON4 a 4 | L341, | 
Gaba | 1 (gS3)+2 56552) | 
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substituting these four values into equation 4.22 


K, (2, v-of fel {re | Lyhg (5 ss at 2b (ao 


b2y 4. C2, C3 
HL, (52) “4, Lyd (24. 1 LGR GP 


Sa. 
Hh DHL, (54) “44, On oT} dxdy 


simplifying again, through the use of Table 4.1 


K, ( 2,4) = bab3te 9¢3+2b1b3+2c4C3tb,b9tc Hope ey 


As can be seen, terms in the K, matrix can be quite lengthy 
and require considerable time to derive. On a more positive 
note though, this square eae is symmetric and thus only 
half the terms need be calculated manually. 

Next, attention will be focused on the Ky and K, 
matrices. These two arrays may be considered together 
since the only difference between the two of them is that 
b values are associated with ka and c values with KE 


Otherwise, they are identical. These two matrices were 


given as 
: P 
K,(i,j) = (ONG NL” )dxd 4.23) 
ot Js ome , 
SL 
and 
Ka(i,j) = y (ght N. )dxdy (4.24) 
‘ eine y 


Taking the (3,6), equation 4.23 becomes 


K,(3,6) = f ea Na) dxdy 


(e) 
where 


a b3 b 
CG .o== if [ Bd he | L,dxdy 


= 2. (L1L,b.+L Speer 
x aS73 oe Tl? aa 


once again, using Table 4.1 


- &£ 
K, (3, 6) Sg De) 


and consequently 


K,(3,6) = (2c, 45) 


Following the same procedure throughout each of these 3x6 


matrices, complete  K, | and | x, | are | 


“P 
by b, +2b, 0 b, bd, 0 bee, 


QO b,+b 


<2 
[K,] =2]o 2d +b, b2 byt2b, 3tb, 
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and also 


0) oe 0) 2cntc., C2 2c, te. 


The next two matrices to be derived, that is | Ko | = 
and | Ks | iz can simply be written down by inspection of 


the two above arrays. Thus 


2 igs Mage |e: 

0 b 0 

- Z 

2 6 

op 

bytb, by +2b, 2b, +b, 
0 0 b, 
b,+2b, bath 2b, +b, 


while its counterpart is then 
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1 
c tc, 2c, te, c, te, 
0) Cy 0 
K oa +c +2 Z 
| 3 ae Co 3 i) C4 e ees 
0 0) C4 
c,t2c., SR 2c, te, 


Finally, the last elemental matrix to be analyzed is the 


one associated with the time-dependent nodal parameters. In 


subsection A this matrix was given as [ cp | . For conven- 
ience here, let [ cp | = IK. | , then 
Ki, 35) — oF ey (4.25) 
ola 


with both i and j running from one to six. Consider, for 


example, point (1,5) 


K,(1,5) = J xrgent 
of 


substituting from equation 4.19 
se DD, =f (2042-14) (22-14) dxdy 
AG Na 


7 EN) «2 
> f (ay L,~-2L,L3.-2L,°L,+.,L,)dxdy 
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then from Table 4.1, using (&%+4+%) four separate times, 


this term reduces rather easily to 


Factoring out a constant of 


— , the total K. matrix takes 
on the form 


Ls] - Le] ats 


16 -4 16 0 a 


Which is also a symmetric matrix, thereby allowing faster . 
derivation of the individual terms with less chance of 


numerical error. 


C. STRUCTURE OF COMPUTER PROGRAMS FOR FLOW ANALYSIS 

A total of three computer programs analyzing two dis- 
tinct test cases of fluid flow problems were employed in 
this thesis. The first was a steady state analysis of 
Couette flow. This involved determining the solution of 
the velocity profiles (linear and nonlinear) in a shear- 


and pressure-induced flow between flat parallel plates. 
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The upper plate slides in the positive x-direction with 

a constant velocity u, while the lower plate remains 
stationary. There is no velocity component normal to the 
plates, that is, v=0 in the y-direction. The governing 


equations for this particular fluid flow are 


Continuity: 

as + st : (4.26) 
Momentum: 

uge + vSu - Jp SP tyv's (4.27) | 

gu + vgve 1 ohiyy'y (4.28) 


For determination of the velocity profile involving only 
linear terms, the left hand side of equations 4.27 and 
4.28 are set equal to zero (no inertia terms). 

A physical representation of the Couette flow analyzed 
is shown in Figure 5. Node and element numbering is the 
same as in Figure 2. A pressure gradient of -3 units is 
directed along the x-axis, i.e., == @3. 

The two remaining programs formed the majority of the 
theoretical portion of this study. They were devised to 
carry out the calculations for the analysis of two-dimensional 


or axisymmetric natural convection heat transfer problems. 
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The first is the lesser-complex steady state approach, 
whereby all transient conditions characteristic of the 
system are assumed to have died out, leaving only those 
steady state parameters remaining to be solved for. Once 
the finite element matrix equations describing the system 
are correctly assembled, a library subroutine (LEQT2F) 
functioning as a linear equation solver is called and the 
desired nodal parameters can be calculated. The second 
program takes into account the previously-neglected time 
dependence of the system by introducing a type of finite 
difference integration scheme to solve the transient por- 
tion of the governing equations. This integration reanmneee 
must be an iterative procedure in order to circumvent the 
problem of nonlinearity similar to that resulting from the 
addition of inertia terms. Furthermore, the integral is 
solved at successive time steps, with time being increased 
until the value of the field variable converges, within 
tolerance, to that of its steady state counterpart. This 
segment of the total equation is then algebraically com- 
bined with the remaining steady state solution to yield 
values of nodal field variables with improved accuracy. 
Gravity acting in the longitudinal direction was taken into 


account for both the two-dimensional and axisymmetric flows. 
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The flow region to be studied is first defined, 
followed by the setting up of coordinate axes (Figure 1). 
The location of the origin of these axes is in most cases 
arbitrary, except that for a problem involving axial 
symmetry the x-axis must coincide with the system axis of 
symmetry. The flow region is then divided into a mesh of 
triangular elements, and the nodal points are numbered in 
the sequence previously described. Once the setting-up 
of the solution demetin is completed, computer analysis of 
the system with its included boundary conditions can be 
initiated. 

Structure of the steady state fluid mechanics problem 
will be discussed in detail, primarily because it comprises 
one entire program and with the addition of the transient 
stiffness matrix elements, accounts for the main program 
in the time-dependent study. 

The program was coded in FORTRAN IV language and begins 
with a series of DIMENSION statements, which set up the 
‘arrays needed in the calculations. As indicated in state- 
ments 0260-0310, storage has been allocated for problems 
with up to 117 nodes; however larger problems can be con- 
- sidered by simply increasing the dimensions of these 


matrices. The limit of problem size is dictated by the 
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core storage of the available machine. Next, the program 
calls for a declaration of the type of program to be 
solved (either two-dimensional or axisymmetric); then the 
appropriate problem label is printed (statements 0360-0410). 
Before proceeding to input the data describing the finite 
element mesh, all the matrix arrays must be initialized 
by setting all terms in these arrays equal to zero. 
Statements 0930-1150 read into the program the node 
numbers and the coordinates of the nodes for the complete 
finite element mesh. Also, the system topology, the ele- 
ment numbers, and the eo of the six nodes associated 
with each element are read. Beginning with statement 1300, 
the velocity, pressure and temperature conditions within 
the solution domain are inserted. Correspondingly, condi- 
tions are specified for the QX, QY, QZC and QZ indices 
where the nodal field variables are unknown. Since the 
solution obtained by the program depends intimately on 
the body of data, the program is queried to print out all 
data that have been input. This enables the programmer 
to check for input data errors. Statement 2970 marks the 
completion of these steps; the program is now ready to 


commence work on a particular fluid mechanics problem. 


66 


The loop to begin calculating the various element 
matrices starts with statement 3010. Once the element 
matrices | ™$ | are computed for one element, they are 
assembled into the master or system stiffness matrix Ea 
by the code followed in statements 7740-7790. The non- 
linear terms appearing in the velocity and temperature 
expressions of the governing equations are formulated in 
statements 4110-5900. An iterative process compares each 
of these terms with a corresponding quantity in the linear 
symmetric [ mg | matrix until they converge in value. It 
is this comparison value that is then assembled with all 
similar values of the other elements to form the system 

[ mm | matrix, which now exhibits or reflects the non- 
linearity. | 

Since each element in the triangular mesh has six (6) 
nodes, the local node numbers are I$ =1, 2, ..096. fhe 
global node numbers for the element are recovered from the 
parameter NODE (K, I$), which was read as input data for the 
element; that is, for element K, the node numbers 
N(1)=NODE(K,1), N(2)=NODE(K,2), etc. were introduced. Then 
the code in statements 7740-/790 loads the terms of the 
elemental matrices into their proper locations in the 


system matrices. Each time that a term of an element 
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matrix is placed in a location in the system matrix where 
another term has already been inserted, this new term is 
added to whatever value is there. A similar loading pro- 
cess takes place for the right-hand-side column vector in 
statements 7810-7940, 

After all the elements have been processed in this 
fashion, the assembled system equations are ready to be 
modified to account for the boundary conditions or 
phenomena. This is done by statements 7980-8060. Thus, 
at the conclusion of statement 8060, the system equations 


possess the form 


Em] {x} - {x18} 


where { Rus } = { xq} = /QXt+gBTm\ , and {x} = 


rity Ss 


Not all of the components of the colum vector { Rus } are 
known because the Q values at nodes where velocity, 
pressure or temperature is specified are unknown; that is, 
at each node i, either u,, Vz, OF T; is known on the one 
side, or Kq,; is known on the other. A similar relation- 
ship exists at the corner nodes for pressure and QZC 


values. The only Q's that can be specifically labeled as 
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heat fluxes are the QZ's, since they directly relate to 
temperature parameters: within the system. 

The only thing that remains to be done now is to call 
a compatible linear equation solver to produce the nodal 
variables sought. In this case, LEQT2F was chosen because 
of its speed and accuracy. 

The following is a list of the symbols and descriptions 


utilized in coding the above program: 


Symbo lL Description 
NCASE interger which specifies the type of 


problem to be solved: 
NCASE=1, 2-D plane problem 
NCASE=2, axisymmetric problem 


NN number of nodes in solution domain 
NNCN | number of corner nodes 

NE number of elements 

XC(1) , YC(TI) | global coordinates of node I 
NODE(J,I) J=1,2,....NE; I=1,2,....6 node 


numbers associated with element J 


NVS(I) node number where velocity or 
temperature is specified 


NPS (TI) node number where pressure is specified 
VELU specified nodal u velocity 

VELV specified nodal v_ velocity 

PNP specified nodal pressure 
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Symbo 1 Description 


INT specified nodal temperature 


NQS (I) node number where a Q value is 
specified; QX and QY are specified 
only at internal nodes, while QZC 
and QZ may be specified at either 
external or internal nodes 


QXNS » . specified nodal value of QX 

QYNS specified nodal value of QY 

QZCNS Specified nodal value of QZC 

QZNS specified nodal value of heat flux QZ 
XC$ (I) , YC$ (I) local coordinates of node I 

TM$ element stiffness matrix 

™ system stiffness matrix 

DEL area of a triangular element 


The program output begins with a statement declaring 
the type of problem to be solved - either nonlinear two- 
dimensional or axisymmetrical. Next, all input data are 
printed and labeled for easy identification. To ensure 
the validity of the solution, the printed input data should 
be carefully checked against the intended input. A state- 
ment following these printed data identifies which nodal 
parameters are associated with which system nodes (re- 
membering specifications of the finite element analysis 


called for a value of both velocities along with a temperature 
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at each node, while a pressure value could be defined only 
at the corner nodes). The complete continuum solution 
follows in the form of a numbered list, in which the 
integer appearing at the far left of this list designates 
the node number, or multiple of it in cases above I35, 
and the figure on the right representing the value of the 


nodal variable in double precision. 


D. NUMERICAL RESULTS 

Complete numerical listings of the field variables for 
both the Couette flow problem and the steady state heat 
transfer problem are shown in the two computer program 
outputs. | 

The velocity profile for the linear Couette flow, i.e. 
node numbers 1-5, 6-10, etc., revealed that the finite- 
element method of analysis agreed with the exact solution 
of this shear-type flow out to the sixth decimal place. 
This is evident by the fact that all five of the FEM points 
lie exactly on the smooth curve depicting the exact solu- 
tion in Figure 15. 

The approximate steady state isotherms of Figure 16 
are directly related to the nonlinear temperatures (node 


numbers 83-117 of the second set of nodal variables) in 


gl 


the heat transfer problem. These isotherms, or constant 
nondimensional temperature lines, vary in value from +1.0 

on the hot temperature wall, to -1.0 on the cold temperature 
wall. The equation used for deriving these values at all 
thirty-five nodal points within the solution domain was 


T-Ty 


eat 
Ty-Ty 


(4.29) 


where T is the nodal temperature and Ty is the mean temper- 
ature of the fluid defined at Ty - a. 

The general shape and relative location of the various 
isotherms within the rectangular enclosure are somewhat 
Similar to those of comparable heat transfer flows involving 
different Prandtl numbers, Grashof numbers and L/D ratios. 
However, due to the relatively low Grashof number of the 
present system (946.4), there is a total lack of a plateau 
in the center region of Figure 16 and the closely packed 
boundary layer flow near the walls is also missing. This 
boundary layer type flow is characteristic of much higher 
Grashof numbers such as those found in the comparative 
examples in | 3 | ; 7] and [ 15] where the Gr; ranged from 
5000 up to 18000. 

Based solely on the thirty-five nodal point temperatures 


available from the solution, the isotherms were sketched as 
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shown in Figure 16. Lacking additional information, the 
shape of the contour lines between such nodes were linearly 
approximated, to a large extent, without speculating as to 
their exact curvature. 

The actual height and width of the enclosure was 
normalized to y* and x*, respectively, for easier inter- 


pretation of the figure. 
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V. EXPERIMENTAL PROCEDURE 


A. ARRANGEMENT OF TEST APPARATUS 

The experimental apparatus was arranged so as to allow 
the study of an essentially two-dimensional fluid flow. 

The major components of the apparatus consisted of the 
test platform, which housed the rectangular enclosure 
(Figure 6), a control system made up of two water circula- 
tors that maintained the vertical copper walls of the test 
platform at desired temperatures (Figure 7), and a large 
(250 mm DIA) plano-convex glass lens for reducing the 
object (8.5 x 1.875 inch vertical rectangular cavity) down 
to a smaller image size that could be completely captured 
on the 4x5 inch holographic plate (Figure 11). 

The rectangular enclosure holding the fluid under 
- investigation was 8.5 inches high, 7 inches long, and 
1.875 inches wide. This test cavity was sandwiched between 
two plexiglas water reservoirs providing constant circula- 
tion by means of manifold connections on their tops and 
bottoms. Hot and cold water drains were located on top 
of the left and right reservoirs, respectively. Similarly, 


on the bottom were the hot and cold water inputs. The 
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inner walls of these two reservoirs were formed by quarter 
inch thick oxygen-free copper plates. Also, these same 
copper plates comprised the principle walls of the rectan- 
gular enclosure, with the ''side walls" being made of plate 
glass, in order to allow visual observations. Six thermo- 
couples were attached to each copper wall and then connected 
to a multichannel recorder for temperature monitoring pur- 
poses. The time needed for each copper plate to reach its. 
respective equilibrium temperature once the water circulators 
had been turned on was approximately 39.8 seconds. For 
comparison purposes, it took the system just under one 
hour (58.1 minutes) for the 50-HB-3520 lubricating fluid 
to attain a steady equilibrium tanpenaime under the same 
experimental conditions. 

An important ——o- imposed by interferometry is 
that the total distance traveled by the object beam must 
be nearly identical with the total path length of the 
reference beam, if the index of refraction throughout is 
uniform. Since the fluid in the rectangular enclosure 
possessed a refraction index of 1.461 and laser light 
along the object beam had to travel through 7 inches of 
this fluid, then the corrected path length through the test 


cavity was 10.23 inches, or a net increase of over 3 inches. 
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With this in mind, the equipment was arranged in a semi- 
elliptical pattern on a heavy table supported at ae 
critical areas by inflated inner tubes. These were to act 
as stabilizing devices. Equipment could not be arranged 
on an exact ellipse due to the fact that a distance of four 
feet, five inches alone was needed from the object to the 
plano-convex lens out of a total table length of eight 
feet. Even so, the turning mirrors were located on the 
apexes of the "shortened" minor axis, the beam splitter at 
one end of the major axis, and the aqueous hologram holder 
at the opposite end (Figure 8). Spatial filters were 
employed to clean up each beam and expand it. A diverging 
lens was inserted just after the object beam spatial filter 
in order to expand this beam to proper size before it 
reached the rectangular test slit. Also, a collimating 
lens was placed between the spatial filter and hologram 
holder along the reference beam. Finally, a large diffusing 
screen made from a piece of developed film mounted on plexi- 
glas and secured in a rigid metal frame, plus the test 
platform itself, were placed es the object mirror and 
the hologram holder (Figures 9 and 10). 

Choosing the correct hologram holder is very important 


in live fringe holography. The reconstructed virtual image 
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must exactly match the original object. Unfortunately, 
after the processing of a hologram, the emulsion on its 
surface tends to dry, thus causing an unwanted displace- 
ment of the virtual image. One procedure that may be used 
to circumvent this problem is to choose a holder that main- 
tains the hologram in aqueous surroundings, such as was 
utilized in this experiment. Also, in order to keep the 
hologram perfectly rigid during and after processing, the 
exposed glass plate was secured in a removable metal frame 
complete with handle. In this way, exact replacement in 
the hologram holder after processing posed no problem. 

Two micrometers built into the holder's top and left side 
were then used for fine adjustment of the hologram. 

The test cavity or rectangular enclosure was filled 
with a very highly viscous fluid (actually a lubricant) 
produced by Union Carbide and known as "UCON" 50-HB-3520. 
This fluid was required to possess physical properties 
such that Rayleigh numbers, based on cavity width, of the 
order of 10-10? could be obtained in the apparatus, at 
accurately measurable temperature differences. Since 
(Ty-T¢) was held constant throughout the experiment, only 
one Rayleigh number was calculated. Its value of 1.0755x104 


was welk within the above tolerance zone. Worth mentioning 


77 


is the generally accepted prediction that above Ra& 2.0x10°, 


the phenomena known as "secondary flows” begin to occur. 
Obviously, such was not the case in this experiment. 

A scribed grid pattern was attached to the back side of 
the test cavity to assist in alignment of the fringes. Two 
water circulators were connected by tubes to manifold nipples 
on the reservoir ends of the test platform. One circulator 
was set to deliver distilled water at 20°C (cold temp.) and 
the other at 25°C (hot temp). By using slide valves, the 
amount of water expended from the circulators could be 
regulated and controlled. 

The experimental procedure was initiated only after the 
entire system had been carefully aligned. The rectangular 
enclosure was allowed to sit undisturbed i 4 period of 
at least several hours to ensure an equilibrium temperature 
state throughout the fluid. Then, a shutter was placed 
directly in front of the beam of a 3 milliwatt, helium- 
neon continuous wave laser serving as the coherent light 
source. After an Agfa-Gevaert Inc. 10E75 holographic 
recording plate was placed in the holder, the shutter was 
opened and the plate was exposed for one and one-half 
seconds (Figure 13). This plate was then removed, developed, 


and replaced in its exact position. Both circulators were 
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» 


then turned on, and a continuous flow of water at 20°C 
and 25°C was allowed to cycle through the reservoirs on 
the test platform. Once fringe lines appear, their 
visibility can be strengthened by following the procedure 
outlined in the next subsection. 

I£ one word could be used to describe the single most 
important factor determining the success or failure of 
this experiment, it would have to be rigidity. All rela- 
tively light-weight gear, such as; the laser, turning 
mirrors, spatial filters, and the plano-convex lens had 
to be weighted down to make them immovable. The test plat- 
form, in which the rectangular cavity was located, was 
sufficiently heavy on its own to preclude it from having to 
be additionally weighted down. The hologram holder already 
came with a very heavy base attached. All connecting de- 
vices, including the tubes transporting heated water to 
the plexiglas reservoirs and plastic sleeves housing the 
thermocouple leads were securely taped together to prevent 
vibration or motion. Any such random vibration would cause 
the fringe patterns to become lost. 

The viewing of these fringes and the subsequent 
collecting of data can be accomplished by positioning the 


appropriate camera in a direct line with the rectangular 
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enclosure, plano-convex lens, and hologram holder (Figure 
11). A television monitor (Figure 12) was employed for 
convenient viewing of the fringe patterns in an area 


adjacent to the experimental set-up. 


B. HOLOGRAPHIC INTERFEROMETRY APPLICATIONS 

Holographic interferometry is an excellent technique 
for developing interference fringe patterns, which may in 
turn be evaluated to quantitatively provide an accurate 
temperature field throughout the domain of the system. 

Heat convection in a rectangular cavity would be diffi- 
cult to analyze empirically. However, by replacing the 
sensors that would ordinarily be used to record temperature 
changes and flow rates with holographic Pee one 
can analyze directly the variation of the density fields 
within the rectangular enclosure. This technique also 
eliminates the inherent change in temperature and flow 
pattern caused by the physical insertion of the sensors 
into the test fluid. 

Real-time holographic interferometry allows a continuous 
flow of information to be recorded at the precise time any 
changes in the observed fluid occur. Single exposure holo- 


grams are utilized with real-time interferometry. A time 
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sequence can be derived for each different viewing position, 
with the use of a single developed hologram. 

Such an exposure technique consists of recording phase 
and amplitude information from an object, in this case the 
rectangular fluid enclosure, onto a holographic plate. The 
recording is accomplished through the use of a reference 
and an object (scene) beam, originating from a single source 
(Figure 13). After processing, the hologram is accurately 
repositioned in its holder. Illuminating the plate with 
the original reference beam results in the primary (virtual) 
image being projected onto the same area as was the object 
(Figure 14). By focusing the object and virtual image beams 
onto a film or focal plane, and then adjusting the system 
so that the two interfering wavefronts (object wave and 
reconstructed wave) coincide, fringes can be produced. 

The hologram now can be finely adjusted to orient the 
fringes in either a vertical or horizontal reference frame. 
Likewise at this time, the fringe patterms can be made to 
appear more visible by varying the beamsplitter to increase 
the intensity of the reference beam while decreasing that 
of the scene beam. If the original object is changed or 
altered in any fashion by the effect of temperature, motion, 


or pressure, an exact superposition will create a reinforcement 
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or cancellation of the intensities of the two waves with 
the result being the establishment of a fringe pattern. 

A dark fringe is produced whenever the difference between 
the object and reconstructed wavefront involves an odd 
factor of 77/2. Bright fringes occur when this difference 
equals an integer value of 277. 

By inserting a camera in-line with the scene (object) 
beam, but on the back side of the holographic plate, one — 
can observe and record live fringe data. This technique 
provides a real time analysis of an unsteady system without 
the need for expensive and time consuming sensors and 
calibration. 

After processing has been completed, problems arising 
from live fringe single exposure holography include; dis- 
placement of the virtual image due to drying emulsions on 
the holographic plate, and non-exact replacement of the 
hologram in its holder. If any relative motion whatsoever 
has transpired between pieces of the experimental equipment 
during or after replacement of the hologram, the fringe 
patterns may be destroyed. 

It was this last problem that caused a particularly’ 
detrimental effect on the experimental results of this 


thesis. Somewhere within the system (apparatus arrangement) 
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there existed a source of motion or a piece of gear 
slightly off the horizontal reference plane that completely 
eliminated these fringe formations almost immediately after 
they evolved. The exact source(s) was never totally iso- 
lated, but the possible choices were reduced down to two, 
the water circulators and/or the support stand of the 


hologram holder. 
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VI. CONCLUSIONS 


The finite-element method was incorporated into steady 
State and time dependent computer programs for analyzing 
laminar convective heat transfer between parallel plates. 
Two sample cases were tested utilizing the general steady 
State program. In each case, values of derived field 
variables compared favorably to either an exact solution, 
in the case of the Couette flow problem (Figure 15); or to 
similar theoretical results, in the case of the heat trans- 
fer problem (Figure 16). The exact solution of the velocity 
profile for Couette flow was obtained by programming the 
analytical expression given by equation 5.5 in [11 | ‘ 

After successful steady state results were achieved, a 
second computer program was then designed to take into 
account the previously-neglected transient behavior of the 
system. A major portion of the steady state fluid mechanics 
problem was interfaced with a series of subroutines, wherein 
the time-dependent terms were to be calculated, to yield a 
total solution.to the governing system of equations and 
associated boundary conditions. Time itself became a 


limiting factor in the completion of this second program. 
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A problem associated with the convergence of the field 
variables remains to be resolved. 

In the experimental phase of this thesis, five attempts 
were made to produce live fringe formations through the use 
of holographic interferometry. In only one of these attempts 
was there observed a momentary interference pattern, corre- 
sponding to the temperature gradient, across the test 
section. This observation lasted approximately three (3) 
seconds after the water heaters/circulators were activated. 

The main factor(s) influencing this inability to acquire 
such live fringes, on film, was the necessary exclusion of 
the hologram holder from the recording plane because of 
limited table length and/or the vibrations generated by 
the water circulators used in the experiment. Either of 
these detrimental conditions could have eliminated com- 
pletely the formation of interference fringe patterns. 

Due to considerable time delay in the acquisition of 
some of the experimental apparatus, no further documenta- 
tion of real time holographic interferometry study could 


be made beyond the previously-mentioned five attempts. 
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APPENDIX A 


FIGURES 


Figure 1 
Rectangular Enclosure 


/ Insulated 


Hot Cold 
Wall | Wall : 
T..=25 C Tn=Z0°C 
: : L=8.5 in. 
FLUID 
(50-HB-3520) 
y 


Insulated 


86 


Figure 2 
Y(I) Discretization of Solution Domain 
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Figure 3 
Global Node-Numbering vs. Local Numbering 
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Figure 4 
Area Coordinates for a Triangle 
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| Figure 5 
F.E.M. Analysis of Couette Flow 
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Figure 6 
Test Platform with 
Rectangular Enclosure 
and Water Reservoirs 


Figure 7 
Water Heaters and 
Circulators with 
Connecting Tubes 


Figure 8 
Table Arrangement (top view) 
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Figure 9 
Apparatus Arrangement with Reference Beam 
Oriented on the Right-hand-side of Table 


Figure 10 
Panoramic View of Experimental Layout 
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Figure 11 
Direct Line-up on Object Beam from 
Test Platform to Recording Camera 


Figure 12 
Television Monitor Used for Convenient Viewing of Fringes 
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Figure 13 
Holographic Recording (top view) 
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Figure 14 
Reconstruction and Recording of 
Interference Patterns 
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Figure 16 
Steady State Isotherms 
(Nonlinear Heat Transfer Problem) 
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APPENDIX B 


BRIEF REVIEW ON CALCULUS OF VARIATIONS 


A fundamental problem in differential calculus is 
extremizing (maximizing or minimizing) a function f(x) for 
a range of the independent variable x. The problem in 
variational calculus is also extremization; however, it is 
concerned chiefly with the extremization of a functional. 
A simple functional, in terms of only one independent 


variable, would have the typical form 


x 
1(G) a fF G, b., 9.) dx 
=< 


2 
where @ = Q(x) and G,. = oa Dees = = 
x 


Summarizing, the two branches of calculus are related in 
that both are concerned with an extremum; one deals with 
number spaces while the other deals with function spaces. 
In variational problems a functional which is charac- 
teristic of the problem is first formed in terms of a 
function (or functions). Then variations of this same 
functional are investigated with a view toward extremizing 


the functional. In some cases this approach results in a 
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closed form, exact solution. But more often, the problem 
must be solved by an approximate method. One such method 

is the Rayleigh-Ritz technique. This approach is preferable 
to the direct application of finite difference methodology 
to solve the differential equation with its associated 
boundary conditions, because the functional can often be 
used to assure convergence of the approximate solution. 

A simple example of variational calculus is the problem 
of finding the plane curve joining two points (X15 y,) and 
(Xo; Yo) which has the shortest length. The solution 
sought here is the function y(x) describing the curve of 
shortest length; the corresponding functional is the length 


of the curve given by 


x 
I(y) = § i+(d¥)* dx 
dx 
oan 


Using the method of variation of calculus implies that of 
all the curves 

Y(x) = y(x) + EW(X) 
which pass through the given end points, the shortest one 
y(x) must be selected. The problem thus reduces to finding 


the function y(x) that makes the integral I(y) a minimum. 
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Generally, in order to minimize the integral 


Ma 
NO 


I(y) = F(x, y, y')dx 


* 


1 
where y' -X, the function y(x) must satisfy the boundary 


conditions and the Euler-Lagrange differential equation 


omm d fort — 
Oy _o— * * 


The previous result could be extended to several 
dependent and independent variables. For example, in 


order to minimize the integral 


I(@) - S free ¥, 9, 9,, O,)dxdy 
A 


in which @, and o, are the partial derivatives of @ with 
respect to x and y, respectively, the general function @ 


must satisfy the Euler-Lagrange differential equation 


OF __9 (OF . & @B 4 
6° dx ‘dG, dy Of, 


in addition to the specified boundary conditions. 
In the past two decades, since the advent of high 
speed digital computers, the variational formulation has 


been quite extensively employed in the fields of structural 
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and continuum mechanics. Important variational principles 
such as least work, minimum strain energy, minimum potential 
energy, and Reissner's variational theorem of elasticity 
have been well developed and are documented in standard 
a However, similar variational principles appli- 
cable to fluid mechanic problems have not been as compre- 
hensively developed. Calculus of variations has, until 
only recently, been utilized sparingly in the field of 


fluid mechanics. 
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TRIS TS 4 2-0 NONLINEAR COUETTE FLOW PROBLEM 


IBANO= 26 


NEQ = 82 


CF NOOES= 35 


CF 


‘ NCQo 


ELEMENTS2 12 


NO e 


CF CCRNEF NOOES= 22 


NCe 


NNVELS# 14 


NNCXY2 21 


NNPS# 12 


SUMMARY OF NCOAL CCCROINATES 
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NUUEDS WMERC BA ANU RT ARE Se ouLiricu 
I NCDE Qx 


f 2) 
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NCCES WHERE FRESSURE IS SPECIFIED 
I NOOE PRESSURE 


NKODOdDVNOHMNSUN 
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@02080e080080808 8 
folelelelololal elelelele) 
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NOCAL VARIABLE IS TRE U-VELOCITY AT NODES t — 35; 
THE V-VELCCITY AT NODES 36 —- 70; 
ANC THE PRESSURES AT NODES 71 -— 826 


TRE FIRST SSCUENCE OF 82 NODAL VARIABLES 
REPRESENTS A LINEAR, STEADY STATE SYSTEM; 

WHILE THE SECOND SET QF THE 82 VALUES CORRESPONDS 
TO A NONLINEAR ANALYSTS QF TRE SOLUTION DOMAIN. 


A PRESSURE GRADIENT IN THE HORIZONTAL SHEAR 
DIRECTION QF ~-3 UNITS IN MAGNITUDE HAS BEEN ADDED 
TO PROOUCE A CURVED VELOCITY PROFILE. 
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NCDE NC. NODE VARIABLES 


1 0210000006000 01 
2 0.1031249SSED 01 
3 0287500000670 00 
4 0253124595610 00 
5 0.0 
a 0210000600000 01 
7 0.10312499990 01 
g 0287500000220 00 
9 053124955510 00 

10 0.0 ' 

11 0.10000000000 01 

12 0210312495960 01 

13 068745S9SS5S8) 00 

14 0. 53124995550 00 

15 0.0 

16 0210000000000 01 

17 0210312500000 01 

18 0.87455555580 00 

1S 025312499S9SD 00 

20 0.0 

21 0.1000C000000 01 

22 0.10312500020 Ol 

23 087500000140 00 

24 0253125006200 00 

25 0.0 

26 0210000000000 01 

27 0.10312506060 01 

28 008745595S86D 00 

29 0.53125000560 00 

30 0.0 

31 0.10000000000 01 

32 0210312500140 01 

33 0. 87455SS887D 00 

34 0. 53125001350 00 

35 0.0 

36 0.0 

37 0.60994265640-15 

38 -0.712£187€18D-42 

39 -0.60954265€40-15 

40 0.0 

41 0.0 


42 0-838102832380-15 
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“3 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
x 
70 
71 
72 
73 
74 
75 
76 
77 
78 
79 
80 
81 
82 


NCOE NOe 


l 


“Ve5ssn soeocoyu-4Z 
-0.88102832380-15 
0.0 
0.0 ; 
0.19352502£2D-18 
-0.11149020410-41 
-0.19252502£20-18 
0.0 
0.0 
0.47096546520-18 
-0-15718222820-41 
-0.4709654£S20-18 
0.0 
0.0 
0.11670534580-17 
-0.218€8403170-41 
-0.11670534580-17 
0.0 
0.0 
0. 2900518 2730-17 
00 2671£796€570-41 
-0.2900£182730-17 
0.0 
0.0 
0.20080511240-17 
-0.3356£885460-41 
-0.20080£11240-17 
0.0 
0.4000000CC0D O01 
0. 40000000000 01 
0.40000000000 ol 
0.30000000000 01 
0.30000000000 01 
0. 30000000000 01 
0.2000C000000 O01 
0.2000C000000 01 
0.2000C000000 01 
0.1000000000D0 01 
0.10000000000 01 
0.10000000000 01 


NOOE VARIABLES 


0.10000000009 Ol 
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iw 


24 


28 


26 


Ve hVUCGOGI“ssu Ub 
0.868958445SD 00 
0.52886443S4D 00 
0.0 | 
010000000000 01 
010251188120 01 
0.8708250£23D 00 
0.5291035280SD 0c 
0.0 
0.10000000000 01 
Oo 10292004470 O1 
0.87078$3553D 00 
0252928248450 00 
0.0 
0.1000000G00D 01 
0.1029££€2570D 01 
0287181955720 00 
0.5295482£46D 00 
0.0 
@.1000000C00D 01 
0. 10297264840 O1 
0.871484366SD 00 
0252970566770 00 
0.0 
G.10000C000000 01 
0.1030036890D 01 
0.8728S57E61D 00 
0252002332500 00 
0.0 
0.1000000C00D 01 
0. 10301734920 Ol 
0.873415902€D 00 
0.5301605230D 00 
0.0 
0.0 
0.94092225430-19 
0. 24063440910-15 
~025690421£$50-19 
0.0 
0.0 
0.12824932130-18 
0. 26110722880-15 
0. 86526782260-15 
0.0 
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4c 
47 
48 
49 
50 
51 
52 
a2 
54 
55 
5&8 
57 
58 
593 
60 
61 
62 
63 
64 
ye! 
66 
67 
68 
6$ 
70 
71 
72 
73 
14 
75 
16 
77 
78 
719 
80 
81 
82 


VeW. 
0. 25283£49350-18 
0. 32006190050-15 
-0.1998€922240-18 
0.0 
0.0 
0.561 2465$840-18 
0. 36351135020-15 
=0.49115243710-13 
0.0 3 
0.0 
0.12815175680-17 
Oe 37806257670-15 
-0612007464650-17 
0.0 
0.0 
0.29622390S2D-17 
0. 32117633760-15 
~0629108457190-17 
0.0 
0.0 
0.20082897E4D-17 
0.15102712170-15 
~0.19508894240-17 
0.0 
0.4000C000000 01 
0.4000000C00D 0o1 
0.4000C000000 01 
0.3000000000D 01 
0.30000000000 01 
0. 30000000000 01 
0.2000C00000D 01 
0.20000000000 01 
0.20000000000 01 
0.10000000000 01 
0.1000C000000 02 
0.1000C000000 0o1 
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STEADY STATE FLUID MECHANICS PROBLEM 


(HEAT TRANSFER) 
THIS IS A 2-D NONLINEAR PROBLEM 


IBAND# 26 
NEC=117 


NGe OF NOOES= 35 

NOe OF ELEMENTS>= 12 

NCe GF CCRNER NODES= 12 
NNVELS= 20 

NNCXY= 15 

NNPS2 6 

NNTS= 10 

NNCZC= 6 


NNCQZ= 25 


SUMMARY OF NODAL CIJCROINATES 


I X(T) Y(T) 
1 0.0 82500 
3 0.0 42250 
3 0.90 0.0 
Ll 0.625 8-500 
13 02625 40250 
15 00625 0.0 
21 16250 82500 
Zo 12250 4e250 
25 1250 0.9 
EN 1.875 8e500 
22 1.975 42250 
35 1.875 0-0 


LISTING OF SYSTEM TOPCLCGY 


ELEMENT NUMBER NODE NUMBERS 
1 1 qt i3 12 
2 1 r4 3 8 
2 3 4 5 9 
4 5 10 LS 14 
5 ll 17 23 22 
6 1) 12 13 18 
7 13 14 is 19 
8 15 20 25 24 
5 21 27 33 32 

10 21 22 23 23 
ll 23 24 25 29 
12 25 30 35 34 


NODES WHERE VELOCITIES ARE SPECIFIED 
I NODE J VELOCITY V VELOCITY 
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NCDES WHERE TEMPERATURE IS SPECIFIED 
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NCQDES WHERE HEAT FLLX CZ IS SPECIFIED 
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= 


NODE HEAT FLUX 


UWS WE OO DAO SU fPor O OOO SW 
ODODANNEWN FS OWAISUPWNFOS OW IO 
SOSCSOSSCSSSOOSSOSSOGO GAGS 
DOOOOOOLOOOOOCOOOOOOOOOCO 


NPD NDI NE er Pr peer rw 
WMI RN PDIDID Fr Pr Pr 


NODAL VARTABLE IS TFE U-VELOCITY AT NODES 1 = 353 
THE V=VELCCITY AT NODES 36 —- 703 

THE PRESSURE AT NODES 71 - 823 

ANO THE TEMPERATURE AT NOOES 83 - 117. 


THE SPECIFIEC WALL PRESSURES 

ARE NORMALIZED TO CNE (1) ATMOSPHERE, 
THAT IS,» 1014000 OYNES/SQeCM 

(ALL PARAMETER VALUES ARE IN CGS UNITS) 


TRE FIRST SEQUENCE CF 117 NODAL VARIABLES 
REPRESENTS A LINEAR,STEACY STATE SYSTEM; 

WHILE THE SECOND SET OF THE 117 VALUES CORRESPONDS 
TG A NGNLINEAR ANALYSTS OF THE SOLUTION OOMAIN. 


NOCE NGe NODE VARIABLE 
l 0.0 
2 0.0 
3 G.0 
4 0.0 
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VeV 

0.0 
-0.15945€01720 
-0. 78826727870 
0015945801720 

0.0 

0.0 
-~0.23649507488D 
-0.120 28962780 
~0. 36455074880 
0.0 

0.0 
-0.12166851530. 
-0.92223087&5D 
~0.121 66891520 

0.0 

0.0 
-0.160518€£15D 
-Qe 10709975830 
-0. 16051865150 

0.0 

9.0 
-0. 1686€13£62D 
—~0e 6458410 240D 
-0. 1685613 £620 

0.0 

0.0 

0.0 

0.0 

0.9 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0. 3287€138200 


01 
01 
01 


01 
090 


01 
02 
O01 


v1 
01 
C1 


00 


~C.86715052020-25 


~0.3287813E2C0 
0.0 
0.0 
0248418167430 


-0.12126888&260- 


00 


00 
22 


igi aa 


be 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
6Q 
6l 
62 
63 
64% 
65 
66 
67 
68 
6S 
70 
71 
UE 
ue) 
14 
75 
76 
7 
78 
3 
80 
81 
&2 
e3 
84 
a5 
86 
&7 
88 
89 
$0 
ak 
92 


“Ue tO4L04014SU UY 
0.0 

0.0 

0. 28532060500 00 
Oo 4210519668D-23 
~022853206050D 00 
0.0 

0.0 

O- 158S2427070-91 
009241557484D-2 3 
-0.1585242707D-01 
0.9 

0.0 
-0.434E077ES70-01 
0.1120825543D-22 
02434E077ES70-01 
0.0 

0.0 

0.9 

0-0 

0.0 

0.0 
0-10140C00C0OD 67 
0-1014CCOOCOD 07 
0-10140C00C00D 07 
02101403 7CE4D 07 
"021014225€57D O07 
0.101403 70E4D 07 
0210139802130 07 
0210128750670 O7 
0210125802150 07 
0210140C0C00D 07 
0-10140C00000 07 
0-10140C00000D 07 
0.25000C00G0D 02 
Q0-2500C0C00C0D 02 
0-2500C00CCOD O02 
0-25000COCCOD O02 
0e25090CO0C00D 02 
Oo241E6EEE67TD O2 
0224166666670 02 
002416666EE7D G2 
0224166666670 02 
Oo241E6E6EETD O2 
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NODE NCe 


onceoiiwu»rtstS w NY ~= 


Vacoodneuaosrososcouw UZ 
023322232230 02 
023322232230 02 
O. 23322232230 02 
0023332232230 02 
0-22500C00C0D 02 
0.22500C00CGD 02 
0.2250C00C00D 02 
0.2250000C00D 02 
0.225C0000CCD 02 
Oe 21666666E7D 02 
Oo21LE66E6EE7D O2 
0. 2166666E67D 02 
0221666666670 02 
Oo 21666666670 02 
Oe 20832332330 02 
0. 208322232230 02 
0. 20832233330 02 
0220833332330 02 
0. 20822332330 02 
0.2000000000D 02 
0. 200GCG00COD 02 
0.2000000000D 02 
0.2000G000000 02 
0220000000000 02 


NODE VARIABLE 


0.90 

0.0 

0.90 

0.0 

0.0 

0.0 
“O.1671SE72270 Ol 
~0081462882870 Ol 
-0.1667709024D Ol 
0-0 

0.0 
-0.3529066201D 00 
~0.1267086543D 02 
~0.3%E47S6S68D 00 
0.0 

0.0 
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af 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 


—“VUVeFSSH (SBN 1ID3UTVE 
-0. 96623861640 O1 
00 88 746546530-01 

0.9 

0.0 
-0.16437S4£500 O1 
-0011188367€10 02 
00 16392742270 Ol 

0.0 

0-0 
—001666266623D 01 
-02068680556020 01 
-0.1659751782D 01 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0239280002440 00 

0.28717151270-03 
—00391E7605S1D 90 

0.0 

0.0 

0259293605440 00 

0.61262434900-0 3 
-0.59071952840 00 

0.0 

0.0 

0237670757430 00 

0. 99603754C60-0 3 
-0. 37314213780 OC 

0.0 

0-0 

0. 63827156120-01 

0. 137604$5030-02 
-0.5987£Z22856D-O1 

0.0 
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“Oa “geo” 


62 -0024295450740-01 
63 001210812751D-02 
64 0.2677262€200-01 
65 0.0 
66 0.0 
67 0.0 
68 0.0 
69 0.0 
70 020 
‘71 0.1614000000D 07 
72 0210140000000 97 
73 0210140000000 07 
14 0.1014038420D G7 
75 0210142616200 07 
76 0.1014038452D G7 
77 0.1012S9S77ED 07 
78 0210138840200 07 
79 021014060754D 07 
20 0210140000000 07 
81 021014G000G0D 07 
82 0.1014000000D 07 
a3 0.25000CCQCOD 02 
a4 0. 25000C00COD 02 
85 0.25000000C0D 02 
£6 0.25000000C0D 92 
a7 0. 25000000000 02 
ag 062411450471D 02 
a9 0024271874287D 02 
90 0024636782520 02 
91 062423258C19D 02 
92 002364452284D 62 
93 Qe 2323056484D 02 
94 0. 23505664810 02 
$5 002423164SE2D 02 
$6 00 23428292€8D 02 
97 0.2227636€50D 02 
98 0.2241£€2105D 02 
95 0022701506210 02 

100 0. 2335£48222D 02 
101 0.2262750S1SD 02 
102 06 2151248262D 02 
103 0621554210590 02 
104 0.2185937501D 02 
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THE 


THE 


THE 


THE 


THE 


THE 


Ve de4OsyusL IU ve © 
Oe 21795465440 02 
002070€230220 902 
0-208060E577D 02 
0220955082280 02 
0021294435330 02 
Ce20918292520 02 
0. 20370258220 02 
0- 20000000000 02 
0e2000000CCOD 02 
0.-20000C060C0D 02 
0-2000000CC0D 02 
0-2000C00CC0D O02 


KINEMATIC VISCCSITY GF FLUID 50-HB-3520 AT 20 DEGREES Ce. 2 
10.90 SCeCM/SEC 


DENSITY OF 50-HB-3520 AT 20 DEGREES Ce = 100595 GM/CC 


CCEFFe OF THERMAL EXPANSION OF 50-HB-3520 AT 20 DEGREES Ce = 
0-002278/DEGREE Ceo 


THERMAL DIFFUSIVITY OF 50-HB-3520 AT 20 OCEGREES Ce = 
0.00102 SQeCM/SEC 


GRASHCF NUMBER (GR(L)) =(G¥B*L**3*(TH=-TC))/V¥*¥2 = 94664 


U VELOCITY FORCING FUNCTION, G*¥B*T (INITIAL), = 652.962 CM/SQS 
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STEADY STATE COUETTE FLOW PRCBLEM 


BAO PION O lM DOOQAXHNASNOM OHO ANM PND OKDORKHCUM SO 
oleleolale(sleolel el fer rie Leal tole lararatiaiatlaialatlatlasradlaslatlasvaslaal 
elololololololololellololofolololololmlolalol~lololelolalolololololole) 
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MODIFICATION OF RHS FOR TM BOUNDARY CONOITIONS 


x< 
sas? | 
os 
x< 
+ 
x< 
~~ 
ee 
= 
Qa. cat 
> = 
eee amt 
=O | 
=—- —Oo0o= 
=Z =—QDO ex 


ce ww @ Oet™ 
sao NOOO Il — 
| iho be OG | Da | 4 
ond) wee OY ome DK Ht LU 
Y) It <me 9 -m 
AW“ eo oo” lO 
teed eee OOK mt et I 
OMN2ZeeS Oe Ot | 
WNww ew Y et tt HO 
OOxKxLTESZSZzTO HOMO 
‘onl on ee Te 4 od ® A OP -e  ) 


w\ 
=< 


na] 


elelelelelelelelelelelelelelelelélel=-\elele 
SQUOAPTINOP-DROSBAIMNS NOH ODKOeEM 
SsfsetesstsataQ\au\11qayuQ\aq qanlam mac 
WDDDDDODDDODDDODD ODODVDVDDWODOOWDOO 
NNNMNNNINNNANNNANNMN NN. Nsw) 
ANNMNMNMNNNNNN NN NN NM. NIM) 
Se rN ed raed eed weed cael ed ced ceed cd cd cad ced eed ced ead od ot 2 Sd 
LLU EL UL UL ER UL UL LL ULL UL UL LLU LLL LL 


OOQOOCOO99 O00 0000000°0 
FADO DHOANUMOLSLNOMOWHOaAN 
OOO OQ OMS FPPSrrspsrsrnnwn 
WMDODDDAD DDD DDDDDDODO®DO 
MNMNNMNUA”NNNNMNNNWNNNNMY) 
NMNMNMNNANNANNNNVIMNNNNN WY) 
a ee ee ee ee ee ee ee ee ee 
WLU ULL UR LER ULL LLL Lk 


Ooo0do 
PTPUVOM 
LOVEDV LON LON LN 
DODACD 
NMMMY 
NMMNWY) 
—_ I! Jj — 
WU Li LL 


— 


OF ELEMENTS=",1I13597/7/3 


Y'gly 5X, 


COORDINATES 9/7/4 
//) 


3 
3 
9 


—_/ @ 
<T—-  —-§-WUl= ® 
tT ee —~ ome (_) preg om O—WIDN 
Omit HNN em Ow @WLOAM KO SI 
—~ BWW) eR ZHOO Sei ath 
OO WWM enrnn © © & OTK & 
@ eHtMOOMM ewe etl OO aX ENxX< 
— OURO HH OMMMOKaZK AMG 
etN etme emi MOLL aeHOwYe @& 
(lL eo Roe © © edt of COWL 
er eOlUuUCYWVi tle © & fe &ex<XY © -QO0M 
OO eOHOOWI> tN tte hme <KXOOMm 
SsHOesen SZWKNYNS w-w SI YOZS & 
Mme ex >We OS HNe Wete = X< 
eel © eZOOZZZZZD~w~ & «DO & & ot 
KW UWIWUWL Le PmP ENS SZ ZO ZZ ZZ ON XM KETO KK 
Zee eRe eR Rem Oct Odd <t<de o e 2 © © @©® © @® MND MELOAN 
ete ee ST IN Ott Ort © OO & oe & oll Oe > a a 
mo Hao | mer SO Ke KR KC KE KK OK OK OK KOKO ON OX 
OSM CW WeYe PL 0S SSS SSS MOP OOO OLD LOAULIIALNIN MI Stl Ne 
uw Qeetwl>D> Tl ee I DD DS OO SS qwew Faerw ewe 
“io thea seZOCN 2Zawwwwwww rR eRe ee eee eee Or Ree ree eR UeF oH 
WICILL ZOHO ff mW CLUNW  t oi di de didi ed didi fd teint EI 
SIR MOU RR eee eR eRe rrr ESSE SSE SSS te P2222 Se° TTWETOS 
ey ee 1 oo peg a pd TF er OOOO oO KOKI oe 
ALOOAWOOCO~HLOL A LLALACLCOONGOOOOVOOOXOOGO0O0OxXGowoono 
QOSORWMOUO Xr SOM ESSSSS SU UU UU UU a Ue Oe Ue ee ee 
om mad _ 


eG Poaite, 


| ee 


~~ 
= 


ND o Mae Sea re eames 


Onn enn eee 
FHHNOmnoOwne 

NANAIMO SF Pin 
aelelelelelol a) o) 
OOMNIAIANAAIACIC 


nn & & & & 


322% 3242324 
8X,*STEADY STATE FLUID MECHANICS PROBLEM* »////) 


2110,2F10.0) 


ITE,2005) 1I,X(T) 


NS © DONVOWMOOMNOTN OMDNOAt NO WO wv 
NN Oo DODO SAAANMMAM MMMOAMOTHP LTO WO O 
‘op Tas) Om NDOOOOQOCOOOCOSO OOQ0O00o 8oO oo O 

Sssetesqetestetes testes ses ett 


Nes at 


eolelelelelelelelelelele (elelelelelele) ooo © oogoo 
KHOANMOALPNOMONES AN CO PUY OP noe ‘) UNV .~OP- 0c et = 
WMODOODOOO0OOOOrP fr eh P= Pe fc © OOO@M Oo 
WMOADDDDDMDDADAOWC WMWDBDODODOO foertelee) © WOOWH fe elee) 
MNMNNNNNINNNNWN NNMNNNNWN MNWMIWY Y NINN NY) 
ANNNMNNMNN NN MN NNNNNNYN MY) Y NM”MNY) NG 
Sao) a! J J.) J.J.) ad ae Se es — a ad aad 
LLhiLihibthibthibtihiber th LLU LiL ib it ot LL iULiuL LL LL UL Lh tL tL 
a Or aA ou 3 
om © Ww Lu et it 
uw N oOo - fad > 
wy e <« 5 ~ ~ 
od o~ tL —~ er 
a f -~ <{ © (am) -— _J 
~ ~ H# WN WO << 
~ ot os Co) Ww oO > ot 
= -— oe“ N e © MN [ tk 
° —Y mA Mm Ta) QO | rc - 
> 1! tl Oe aa fo) - Sac 
og Q + Ze ~ t ve <{ — oe 
© m Oe ~ co if ~~ ~~ ~ 
~ =z I az ~ a mu © So ml pf ~ 
em @ . amr C/O o { ce wh La | % % ~ 
A aed ee OO ed “NO tlie = ©S ea] Ww 3b c ~ 
eo On LL = <q ~<a OVS =. uy se tb ag) —_—I + ~ 
ax wml emis, ~N> eo QO lw ao & t *% © ~ 
ws » Nete LU™ >. ~ 2©& CO —JI= © wm cO e = 
—t ONLI & tt. mj eet) - On - So 2 i t 2&2 6 
ie il! emsaWe. _ eaiz Zz ~ ON D “! OC f =» oO _ — 
—e me ©) LUX © eOuno no we cw i= So "= ee WY) 
Oe UWWUI-F-OD = C€ ™Ok-Liue - Ae tw —_—- Ww iw eke Ne & 
Wo aey~aAdqdNIeN AZIDE © = riz Ut m= 
a. ODMNa WON Ue &@t KID W ow WL <{ te - 2ae Zz 
Y wy WM Le De MOIS J = Oo oO a rs) - D DW D 
ANNA) & { oo [=~ >0O DH ALD 9 4 aI ul Mt 
ti NWS f€£lwo DOr ett YW <t FO > wu ud > on NA W 
w O WIN Uden ~e OM FM — Wor fF ra) ~ mw © WM © 
< MOAWeFGOTOGD WwW eOlqgqei ££ TQ om [ ad —- OF 2 CA oO 
me Me eae tie LTInN Ores It xt = /) MO « > ~w et Oso 
~ eD eX eQ ee ee OWN YN Fe > §-~X OC as) = WO a 
We >< = << DD KY Oo Ws tt OM<f WwW | oc ”Y Cc Ct jim 
* EFOAQINUNSNY WW YNOMIOO>TuL thi mtli WW Oo tu => wo ._.J e 
Ox DAC aNet —- OZAareO A KZ = Ww =I LL Oo Ww {OW 
ZH Well eae ps Nat fuict -—O8 <— > — LL. = SH we 
= WwetwtwWowe CO Fe OkeOMMe O22 & Oo w (am) z_ eA MN 
<—W COWOWONOAMOMOMNKSUY Nee io = Qo = me LILI) 
GO aodrKrOTOOGIO exZ2U NN HKrNX MZ -&— > ee ee) ee Oo bode 
Om 2 2 2S ZN AU Mi O <f - e ~t NO OO LOYD 
tuawmilie tlie tle the FMR OMRIndIuL em IO FS we WL ez “NT were OC 
MH MA aM oO eo tLe tZDuMi Za we NN | uw “WO WL © & wy Wewoo <a 
eu etl o< LL <LI O< tL < o> meZOtK ae ON 2S NZ tl ew ect S>Tli-o> 
Carter rrrrrcT heh eRe Tee Z NOO- ww ©O OF Wa “+o, O 
S tS oS oT oF AOQOYIKACLIOAGwsS Ze gt YY © OM YO = wo DnnO0m 
mai © © © © wtcfe tw Oda ems ro) uw ”) NUP 
Y YemMemAmNier Oto Wiel LLL ellie Ww Wil tui Ut NM UU LWIA Ape 
QO COO OD OD 0D om ZIIDUIE Ye MOOF OE SN EE OREN ORK OF Oe eR 
Oo ws9Dx<Ox<XOx<XO™XM eo UWNkKe Mumm we YNNe NZe ce etlle Ge ee dt <c 
Z wZOZNZHA™NMK oD eZ? TZ wo KM ee et) & were Oe Ce Us) eS” a 
es Ne eo ewe eo &Mode i wlU>x<iurOxX< Moc DN 8 ex xX lS NK mex AH<_t 
So Fen en ener ane ciInwn ZMx< st-NOe GUO SO TO ONO eNOe0. 


> HM oC NOC NOC NOC NM) oe 


Ak etlitli —ilaMmAe Ned HN ef — 


Oma HG eZ 


W><e A I UW oo ft NYKO ot © oN ew aA e OQ ee YN Ff LI 
eAlee © ee eae ee Two en efter @eOrN™N CNN ee SNC 

~ erm NN NN MN ERTZNWTO KN KC aNNOA™N eiNOnrNON LNT NI 

~~ LAN SOR MND Fe IN EFRON DOMRANN ene NWO™ ©O™M @ON TE ~ 

wu eeers e  O ®  rwrw ww w s O~ ~~ LL ewe ows Om we ew oo 

RR Re eR RR ®t ne oe REE EE UR OF OF HORT EAE © o & 

It tg gf ge ft fein XM eM eM eeae(efitditie Cnet 0 HCOTOd KX 

[STE ETE OSE C ECC TV] MHMNWONTSMNOMNMOSDSATSAAS eVPSTUi eZee esl ZVVMWNWNGA 
ew ane wo wo Cen ewww enenw ONO CCA OKCKCOCOKCORMAHA eK & & &~OO 
0 0C 0 OO O DONNNONNSNCOOOCOCOCOWRO CO AO AO CGHNONNNEZ 
Ui Wilh kh hk OULU eo NULL NNN ULEL LEL WLLL OOOAILL OW ONLI AILL SILL OWL NN Sl 
am +S t =e = HAI AAI aN Sf AN AN SS ANIM 

oO NnoOoAtN M WO aan MONS AW So ww © me ww 

~- ~—Oo © ODO © WO oO KHOOANAIN my Om ~ rT WwW WY 
oo0o0o0 8OlCUOFlUCUOO © oleolelelslol eo) So oO © oOo Oo Oo 

SS ss S| = Ft tH aoa et NANIAICIOCIO N N a | oN i 9 


PZ 


TIME-DEPENDENT FLUID MECHANICS PROBLEM 


ANOTIOM DHOHKHAIN SMO DHOHNM TNO DHAOHINMOSTLVO 
ODVOVW OCC OO et Aa AH HAH AMAA AIANANNNANAAUN MM AMM Mam 
OO OSGASOGOQq0c GO Go qoocdqgqcoodWdood0d0odo00 0000 [0 
olelelealelelealelelelelelelealelelalelelelelalelelelealealealeolealealealolalal@) 
CLIALPRRSALO RAS ULILLIRLIRRILRI AREER ELSA E ILLIA LAPUA APART PPR LR LL 
> Si > > Si > > a a a a a a a a 
PN PN a ed ed pet Et i fe ps a i eS st 


teh REE REE RE EEE I- REE EERE EERE ER EEE EEE EE 


EAL*8 (A-H,0-Z,$) 
/5/ 


} 
I 
3 
S 
OF CORNER NODES) 


PW pe ZS 
pm = OC <I DU 
mm Of ZO 
OOMWZZ 

OCOowa+ W 
JI jw zw 
UL OC LU 2 

>>. + <a f— 


za” 
7>0.-- ZWD 

+e 
ULL LL LL Li) 


TrIrrms~is 
a a a 


Bi 3 


( 
( 
»NQIS (117) 


Ne ee 
oad wer ee OO) 
we /) md eer ne 
C) 8 ee om HF WF om “Sf 
eS ort OO etd & oe 
aZnQAINIDS Alm 
we olf © eo om off 
es AA] aA eed Ht = 
QAM At N01 NA et et 
O fae wee td C/) HF 0H HP 
SNOW GTA O~ et 
Wise << ah OK et 
i ~ 
MA Z2SZZZZZZZO 
CLODODOOOOOO\O 
Le a et tt tg 
ZNNNNNNNWYN 2 
222222220 
Tt <I WW Wl 
ed ae ee ae a a i 
I, <X, pnd 4 P= 4 Pt 9 BG 99 =) 
OOQODOOORQOO0W 


(117,132) 
M(LL7,117)sRHS(LL7) 


CLUDING NON-LINEAR TERMS, 


a 


0 

) 

2 
RITE,2020) 


E 
E 
W 
W 
E 


THE FIRST PART OF THE PROGRAM CAN BE CONSIDERED AS AN INPUT ROUTINE 
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